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ABSTRACT 

We measure the anisotropy of dark matter flows on smaU scales 500 kpc) in the near en- 
vironment of haloes using a large set of simulations. We rely on two different approaches to 
quantify the anisotropy of the cosmic infall: we measure the flows at the haloes' virial radius 
while describing the infalling matter via fluxes through a spherical shell; we measure the spa- 
tial and kinematical distributions of satellites and substructures around haloes detected by the 
subclump finder ADAPTAHOP first described in the appendixiBl The two methods are found 
to be in agreement both qualitatively and quantitatively via one and two points statistics. 
The peripheral and advected momentum is correlated with the spin of the embeded halo at 
a level of 30% and 50%. The infall takes place preferentially in the plane perpendicular to 
the direction defined by the halo's spin. We computed the excess of equatorial accretion both 
through rings and via a harmonic expansion of the infall. 

The level of anisotropy of infalling matter is found to be ~ 15%. The substructures have 
their spin orthogonal to their velocity vector in the halo's rest frame at a level of about 5%, 
suggestive of an image of a flow along filamentary structures which provides an explanation 
for the measured anisotropy. Using a 'synthetic' stacked halo, it is shown that the satellites' 
positions and orientations relative to the direction of the halo's spin are not random even in 
projection. The average ellipticity of stacked haloes is 10%, while the alignment excess in 
projection reaches 2%. All measured correlations are fitted by a simple 3 parameters model. 
We conclude that a halo does not see its environment as an isotropic perturbation, investi- 
gate how the anisotropy is propagated inwards using perturbation theory, and discuss briefly 
implications for weak lensing, warps and the thickness of galactic disks. 

Key words: Cosmology: simulations. Galaxies : formation. 



1 INTRODUCTION 

Isotropy is one of the fundamental assumptions in modern cos- 
mology and is widely verified on very large scales, both in large 
galaxies' surveys and in numerical simulations. However on scales 
of a few Mpc, the matter distribution is structured in clusters and 
filaments. The issue of the anisotropy down to galactic and clus- 
ter scales has long been studied, as it is related to the search 
for large scale structuration in the near-environment of galax- 
ies. For example, both observational study (e.g. West ( 1994), 
IPlionis & Basilakost (2002). Kitzbichler & Sauret (2003)) and nu- 
merical investigations (e.g. Fal tenbacher et aL C2002) ) showed that 
galaxies tend to be aligned with their neighbours and support 
the vision of anisotropic mergers along filamentary structures. On 
smaller scales, simulations of rich clusters showed that the shape 
and velocity ellipsoids of haloes tend to be aligned with the distri- 
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button of infalling satellites which is strongly anisotropic iTormenI 
However the point is still moot and recent publications did 
not confirm such an anisotropy using resimulated haloes; they pro- 
posed 20% as a maximum for the anisotropy level of the satellites 
distribution jVitvitska et"d] <2002h ). 

When considering preferential directions within the large 
scale cosmic web, the picture that comes naturally to mind is 
one involving these long filamentary structures linking large clus- 
ters to one other. The flow of haloes within these filaments can 
be responsible for the emergence of preferential directions and 
alignments. Previous publications showed that the distributions 
of spin vectors are not random. For example, haloes in simula- 
tions tend to have their spin pointing orthogonally to the fila- 
ments' direction ( Falteribacher et al. (2002)). Furthermore, down 
to galactic scales, the angular mome ntum remains mainly aligned 
within haloes iBuUock et alJ i200lh ). Combined with the results 
suggesting that haloes' spins are mostly sensitive to recent infall 
jvan Haarlem & van de Wevg^ ilW 3|)), these alignment proper- 
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ties fit well with accretion scenarii along special directions : angular 
momentum can be considered as a good marker to test this picture. 

Most of these previous studies focused on the fact that align- 
ments and preferential directions are consequences of the formation 
process of haloes. However, the effects of such preferential direc- 
tions on the inner properties of galaxies have been less addressed. It 
is widely accepted that the properties of galaxies partly result from 
their interactions with their environments. While the amplitude of 
the interactions is an important parameter, some issues cannot be 
studied without taking into account the spatial extension of these 
interactions. For example, a warp may be generated by the torgue 
imposed by infalling matter on the disk jOstriker & Binnev ( 198^. 
iLopez-Corredoira et alJ i2002h ') : the direction but also the ampli- 
tude of the warp are a direct consequence of the spatial configu- 
ration of the perturbation. Similarly, it is likely that disks' thick- 
ening due to infall is not independent of the incoming direction 
of satellites (e.g. Ouinn et al. ( 1993), VelazQuez & White ( 19^, 
iHuang & Carlb erg ( 1997)). 

Is it possible to observe the small-scale alignment ? In par- 
ticular, weak lensing deals with effects as small as the l evel o f 
detected anisotropy (if not smaller) (e.g. Hatton & Nininl i200lh . 
ICroft & Metzler (2000), Heavens et al. (2000)), hence the impor- 
tance to put quantitative constraints on the existence of alignments 
on small scales. Therefore, the present paper also addresses the is- 
sue of detecting preferential projected orientations on the sky of 
substructures within haloes. 

Our main aim is to provide quantitative measurements to study 
the consequences of the existence of preferential directions on the 
dynamical properties of haloes and galaxies, and on the observa- 
tion of galaxy alignments. Hence our point of view is more galac- 
tocentric (or cluster-centric) than previous studies. We search for 
local alignment properties on scales of a few hundred kpc. Using a 
large sample of low resolution numerical simulations, we aim to ex- 
tract quantitative results from a large number of halo environments. 
We reach a higher level of statistical significance while reducing 
the cosmic variance. We applied two complementary approaches to 
study the anisotropy around haloes: the first one is particulate and 
uses a new substructure detection tool ADAPTAHOP, the other one 
is the spherical galactocentric fluid approach. Using two methods, 
we can assess the self-consistency of our results. 

After a brief description of our set of simulations (§2), we 
describe the galactocentric point of view and study the properties 
of angular momentum and infall anisotropy measured at the virial 
radius (§3). In (§4) we focus on anisotropy in the distribution of 
discrete satellites and substructures and we study the properties of 
the satellites' proper spin , which provides an explanation for the 
detected anisotropy. In (§5) we discuss the level of anisotropy as 
seen in projection on the plane of the sky. We then investigate how 
the anisotropic infall is propagated inwards and discuss the possi- 
ble implications of our results to weak lensing and to the dynam- 
ics of the disk through warp generation and disk thickening (§6). 
Conclusions and prospects follow. The appendix describes the sub- 
structures detection tool ADAPTAHOP together with the relevant 
aspects of one point centered statistics on the sphere. We also for- 
mally derive there the perturbative inward propagation of infalling 
fluxes into a collisionless self gravitating sphere. 



2 SIMULATIONS 

In order to achieve a sufficient sample and ensure a convergence 
of the measurements, we produced a set of ~ 500 simulations. 
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Figure 1. The number density, dN/dM, of haloes with a mass M (within 
dM) detected in our simulations at redshift zero (solid line) compared to 
the distribution expected by Press-Schechter (dashed line). The agreement 
is good for masses between 3 ■ lO^Af© and 3 ■ IQ^^Mq. 

Each of them consists of a 50 /i^^ Mpc'' box containing 128"' par- 
ticles. The mass resolution is 5 ■ lO'^M©. A ACDM cosmogony 
(Vim. = 0.3, r^A = 0.7, h = 0.7 and erg = 0.928) is implemented 
with different initial co nditions. Thes e initia l conditions were pro- 
duced with GRAFIC fSertschi ngeil i200ll) ) where we chose a 
BBKS (Bardeen et al. ( 1986)) transfer function to compute the ini- 
tial power spectrum. The initial conditions were u sed as inputs 
to the parallel version of the treecode GADGET (iSpringel et alJ 
J200ll) ). We set the softening length to 19 kp c. The halo detec- 
tion was performed using the halo finder HOP iEisenstein & Hull 
( 1998)). We employed the density thresholds suggested by the 

authors (Aoutcr = 80, ^saddle ~ 2.55outcr, 5poak = S.^outcr) 

As a check, we computed the halo mass distribution. It is shown 
in Fig. [T] and compared to the Press-Schechter mass function 
JPress feSchechteJ il974l) ). The measured distribution is in agree- 
ment with the theoretical curve up to masses ~ 3 ■ IO^'^A/q (100 
000 particles), which validates our completeness in mass. 

As an other means to check our simulations and to evaluate 
the convergence ensured by our large set of haloes, we computed 
the probability distribution of the spin parameter A', defined as 
( Bullock etal. ( 2001.) ): 

J 



X' 



V2MVR2 



(1) 



Here J is the angular momentum contained in a sphere of virial 
radius i?200 with a mass M and V'^ = GM/R20Q. The measure- 
ment was performed on 100 000 haloes with a mass larger than 
5 • 10^^ A^Q as explained in the next section. The resulting distribu- 
tion for A' is shown in Fig.|2| Th e distribution P(A') is well fitted 
by a log-normal distribution fe.g. lBuUock et alJi200ll) V 



P(A')dA' = 



A'V2^cr 



exp 



_l£!(VZA|0^ dA'. 

2cr2 



(2) 



We found Aq = 0.0347 ± 0.0006 and cr = 0.63 ± 0.02 as 
best-fit values and they are consistent with parameters found by 
IPeiraniet alJj2003h (A' = 0.035 and cr = 0.57) but our value of a 
is slightly larger. However, using cr = 0.57 does not lead to a sig- 
nificantly different result. The value of cr is not strongly constrained 
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Figure 2. The distribution of the spin parameter A' defined as A' = 
J/ {\/2MV R200 ) computed using 100 000 haloes with a mass greater than 
5 • 10^^ Mq. The distribution can be fit with a log-normal function with pa- 
rameters A[, = 0.0347 ± 0.0006 and o- = 0.63 ± 0.02 (red line). The 
curve parametrized by Ap = 0.035 and cr = 0.57 is also shown (black 
line). The two results are almost concurrent, indicating that the value of a 
is not so strongly constrained using a log-normal distribution. 

and no real disagreement exists between our and their best-fit val- 
ues. The halo's spin, on which some of the following investigations 
are based, is computed accurately. 



3 A GALACTOCENTRIC POINT OF VIEW 

The analysis of exchange processes between the haloes and the in- 
tergalactic medium will be carried out using two methods. The first 
one can be described as 'discrete'. The accreted objects are explic- 
itly counted as particles or particle groups. This approach will be 
applied and discussed later in this paper. The other method relies on 
measuring directly relevant quantities on a surface at the interface 
between the halo and the intergalactic medium. In this approach, 
the measured quantities are scalar, vector or tensor fluxes, and we 
assign to them flux densities. The flux density representation al- 
lows us to describe the angular distribution and temporal coherence 
of infalling objects or quantities related to this infall. The formal re- 
lation between a flux density, ^^(fJ), and its associated total flux 
through a region S, $, is: 

^= I vj{n) ■ dn, (3) 

J. 5 

where fl denotes the position on the surface where zu is evaluated 
and d$7 is the surface element normal to this surface. Examples of 
flux densities are mass flux density, pVr, or accreted angular mo- 
mentum, pVrL. In particular, this description in terms of a spherical 
boundary condition is well-suited to study the dynamical stability 
and response of galactic systems. In this section, these fields are 
used as probes of the environment of haloes. 

3.1 Halo analysis 

Once a halo is detected, we study its environment using a galacto- 
centric point of view. The relevant fields 1x7 (f2) are measured on 



the surface of a sphere centered on the halo's centre of mass with 
a R200 radius (where 3M/(47r7?ioo) = 200p) (cf. fig |3j. There 
is no exact, nor unique, definition of the halo's outer boundary and 
our choice of a -R200 (also called the virial radius) is the result of 
a compromise between a large distance to the halo's center and a 
good signal-to-noise ratio in the spherical density fields determina- 
tion. 

We used 40 x 40 regularly sampled maps in spherical angles 
fl = (f)), allowing for an angular resolution of 9 degrees. We 
take into account haloes with a minimum number of 1000 parti- 
cles, which gives a good representation of high density regions on 
the sphere. This minimum corresponds to 5 • IO^^A/q for a halo, 
and allows us to reach a total number of 10 000 haloes at z=2 and 
50 000 haloes at z=0. This range of mass corresponds to a some- 
what high value for a typical L, galaxy but results from our com- 
promise between resolution and sample size. Detailed analysis of 
the effects of resolution is postponed to Aubert & Pichon 1 2004). 

The density, p{ft), on the sphere is computed using the par- 
ticles located in a shell with a radius of R200 and a thickness of 
^200/10 (this is quite similar in spirit to the count in cell tech- 
niques widely used in analyzing the large scale structures, but in 
the context of a sphere the cells are shell segments). Weighting the 
density with quantities such as the radial velocity or the angular 
momentum of each particle contained within the shell, the associ- 
ated spherical fields, pvr(il) or pL(J7), can be calculated for each 
halo. Two examples of spherical maps are given in Fig.|3| They il- 
lustrate a frequently observed discrepancy between the two types 
of spherical fields, p(f2) and pVi.(fl). The spherical density field, 
p(f2), is strongly quadrupolar, which is due to the intersection of 
the halo triaxial 3-dimensional density field by our 2-dimensional 
virtual sphere. By contrast the flux density of matter, pvr(f2), does 
not have such quadrupolar distribution. The contribution of halo 
particles to the net flux density is small compared to the contribu- 
tion of particles coming from the outer intergalactic region. 

3.2 Two-points statistics: advected momentum and halo's 
spin 

The influence of infalling matter on the dynamical state of a galaxy 
depends on whether or not the infall occurs inside or outside the 
galactic plane. If the infalling matter is orbiting in the galactic 
plane, its angular momentum is aligned with the angular momen- 
tum of the disk. Taking the halo's spin as a reference for the direc- 
tion of the 'galactic' plane, we want to quantify the level of align- 
ment of the orbital angular momentum of peripheral structures (i.e. 
as measured on the virial sphere) relative to that spin. The inner 
spin S is calculated using the positions-velocities (rpart, Vpart) of 
the particles inside the -R200 sphere in the centre of mass rest frame 
(ro,vo): 

S = ^ (rpart - ro) X (Vpart - Vq). (4) 
part 

ro is the position of the halo centre of mass, while vo stands for the 
average velocity of the halo's particles. This choice of rest frame 
is not unique; another option would have been to take the most 
bounded particle as a reference. Nevertheless, given the considered 
mass range, no significant alteration of the results is to be expected. 
The total angular momentum, Lt (measured at the virial radius, 
R200) is computed for each halo using the spherical field pL(f2): 

Lt = / pL(n)df2. (5) 

J47r 
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Figure 4. Excess probability, 1+$ls (^)> of the angle, d, between the halo's 
spin (S) and the angular momentum (Ly for total, or La for accreted) 
measured on the virial sphere using the fluid located at the vrrial radius. 
Here represents the total angular momentum measured on the virial 
sphere (corresponding curve : solid line, circle symbols) and the total 
accreted angular momentum measured on the sphere (corresponding curve: 
dashed line, diamond symbols). The error bars represent the 3o" dispersion 
measured on subsamples of 10 000 haloes. The congelation takes into ac- 
count the uncertainty on the angle determination due to the small number 
of particles at the virial radius. Here £,ijsifi) = would be expected for an 
isotropic distribution of angles between S and L while the measured distri- 
butions indicate that the aligned configuration (9 ~ 0) is significantly more 
likely. The two excess probability distibutions are well-fitted by Gaussian 
functions (red curves: see main text). 



(7) 



Figure 3. A galactocentric point of view of the density field, p(f2) (top) 
and of the flux density of mass, pvr(n), surrounding the same halo (bot- 
tom). This measurement was extracted from a ACDM cosmological sim- 
ulation. The considered halo contained about 10^^ Mq or 2000 pailicles. 
The high density zones ai'e darker. The density's spherical field shows a 
strong quadrupolar component with high density zones near the two poles 
while this component is less important for the mass flux density field mea- 
sured on the sphere. This discrepancy between the two spherical fields is 
common and reflects the shape of the halo as discussed in the main text. 



The angle, 9, between the spin of the inner particles S and the total 
orbital momentum Lt of 'peripheral' particles is then easily com- 
puted: 

Lt ■ S , 



LtIISI 



(6) 



Measuring this angle 9 for all the haloes of our simulations al- 
low us to derive a raw probability distribution of angle, dr{9). 
An isotropic distribution corresponds to a non-uniform probability 
density diso{9) . Typically diso is smaller near the poles (i.e. near the 
region of alignment) leading to a larger correction for these angles 
and to larger error bars in these regions (see fig.|4}: this is the con- 
sequence of smaller solid angles in the polar regions (which scales 
like ~ sin 9) than in equatorial regions for a given 9 aperture. The 
true anisotropy is estimated by measuring the ratio: 



Here, 1 + ^ls(^) measures the excess probability of finding S 
and Lt away from each other, while (,ls{9) is the cross correla- 
tion of the angles of S and Lt- Thus having ^ls{9) > (resp. 
^LS (9) < 0) implies an excess (resp. a lack) of configurations with 
a 9 separation relative to an isotropic situation. 

To take into account the error in the determination of 9, each 
count (or Dirac distribution) is replaced with a Gaussian distribu- 
tion and contributes to several bins: 



9o)^J^i9o,ao) = 



exp 



{9__eol 

2al 



(8) 



o"o \/27r 

where M stands for a normalized Gaussian distribution and where 
the angle uncertainty is ap proximated by gp ^ (i n/N)^^^ using 
N particles as suggested bv lHatton & NiniJ i200ll) . If A'^„ is equal 
to the number of particles used to compute pL{fl) on the virial 
sphere and if Nh is the number of particles used to compute the 
halo spin, the error we associated to the angle between the angular 
momentum at the virial sphere and the halo spin is: 



(9) 



ao = V(47r/7V„) + (47r/iV,0 ~ V(47r/iV,), 

because we have A*'„ ^ Nh- Note that this Gaussian correction 
introduces a bias in mass: a large infall event (large A*'„, small ao) 
is weighted more for a given 9a than a small infall (small A'^„ , large 
(To). All the distributions are added to give the final distribution: 



dr{9) = '^JV{9p,ap 



(10) 
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where Np stands for the total number of measurements (i.e. the 
total number of haloes in our set of simulations). The corresponding 
isotropic angle distribution is derived using the same set of errors 
randomly redistributed: 

di,,{e) = ^N{&;°,ap). (11) 

p 

Fig. El shows the excess probability, 1 + ^ls(S), of the angle 
between the total orbital momentum of particles at the virial radius 
Lt and the halo spin S. The solid line is the correlation deduced 
from 40 000 haloes at redshift z = 0. The error bars were deter- 
mined using 50 subsamples of 10 000 haloes extracted from the 
whole set of available data. An average Monte-Carlo correlation 
and a Monte-Carlo dispersion a is extracted. In Fig.|4| the symbols 
stand for the average Monte-Carlo correlation, while the vertical 
error bars stand for the 3(t dispersion. 

The correlation in Fig.|4|shows that all angles are not equiva- 
lent since ^ls(6) 7^ 0. It can be fitted with a Gaussian curve using 
the following parametrization: 



Hi 



27ra3 



exp 



-{e-a2f/{2al)\ +aA. 



(12) 



The best fit parameters are Oi = 2.351 ± 0.006, a2 — —0.178 ± 
0.002, as = 1.343±0.002, 04 = 0.6691±0.0004. Theinaximum 
being located at small angles, the aligned configuration, LtS ~ 0, 
is the most enhanced configuration (relative to an isotropic distri- 
bution of angle 9). The aligned configuration of Lt relative to S 
is 35% (5ls(0) = 0.35) more frequent in our measurements than 
for a random orientation of Ly. As a consequence, matter is pref- 
erentially located in the plane perpendicular to the spin, which is 
hereafter referred to as the 'equatoriar plane. 

The angles, are measured relative to the simulation 

boxes z-axes and x-axes and not relative to the direction of the spin. 
Thus we do no expect artificial Lt-S correlations due to the sam- 
pling procedure. Nevertheless it is expected on geometrical ground 
that the aligned configuration is more likely since the contribution 
of recent infalling dark matter to the halo's spin is important. As a 
check, the same correlation was computed using the total advected 
orbital momentum: 



(13) 



The resulting correlation (see Fig.|4j is similar to the previous one 
but the slope toward small values of 9 is even stronger and for ex- 
ample the excess of aligned configuration reaches the level of 50% 
(^Ls(O) ~ 0.5). The correlation can be fitted following Eg . 11 21 with 
01 = 3.370 ± 0.099, 02 = -0.884 ± 0.037, as = 1.285 ± 0.016 
and a4 = 0.728 ± 0.001. This enhancement confirms the relevance 
of advected momentum for the build-up of the halo's spin, though 
the increase in amplitude is limited to 0.2 for 9 = 0. The halo's in- 
ner spin is dominated by the orbital momentum of infalling clumps 
(given the larger lever arm of these virialised clumps) t hat have just 
passe d through the virial sphere, as suggested by Vitvitska et alj 
12002) (see also appendix ID! . It reflects a temporal coherence of 
the infall of matter and thus of angular momentum, and a geomet- 
rical effect: a fluid clump which is just being accreted can intersect 
the virtual virial sphere, being in part both "inside" and "outside" 
the sphere. Thus it is expected that the halo spin S and the momenta 
Lt and La at the virial radius are correlated since the halo's spin 
is dominantly set by the properties of the angular momentum in its 
outer region. The anisotropy of the two fields Lt and La do not 
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Figure 5. Top: Normalized probability distributions (PDF) of the excess 
of equatorial infall, 5m , measured at the virial radius. The quantity 1 -|- 5m 
stands for the ratio between the flux of matter through the equatorial sub- 
region of the -R200 sphere and the average flux of matter thi'ough the whole 
R200 sphere. The equatorial sub-region is defined as being peipendicular 
to the direction of the halo's spin. It formally corresponds to the top-hat- 
smoothed mass flux density contrast. <5„i = is expected for an isotropic 
infall of matter thi'ough the virial sphere. The average value of 5m is al- 
ways greater than indicating that the infall of matter is, on average, more 
important in the direction orthogonal to the halo's spin vector than in other 
directions. Bottom: the antisymmetric part of the 5m distribution. Being 
positive for positive values of 5m , the antisymmetric part of the 5m distri- 
butions shows that accretion in the equatorial plane is in excess relative to 
the one expected from an isotropic accretion of matter 



have the same implication. The spatial distribution of advected an- 
gular momentum. La, contains stronger dynamical information. In 
particular, the variation of the halo-l-disk's angular momentum is in- 
duced by tidal torques but also by accreted momentum for an open 
system. For example the anisotropy of LpVr should be reflected in 
the statistical properties of warped disks as discussed later in sec- 
tions 3rid|6!23 
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Figure 6. The redshift evolution of (<5m). The average (■) is performed on a 
set of 40 000 haloes at z = and 10 500 haloes at z = 1.8. The error bars 
stand for the error on the estimation of {Sm ) with A = cr {8m ) / VN. N is 
the number of haloes needed to compute {5m, ) ■ {Sm ) is always positive and 
indicates an excess of accretion in the equatorial plane. This redshift evo- 
lution can be fitted as {5m)(z) = 0.0161(±0.0103)z + 0.147{±0.005). 
This excess is detected for every redshift smaller than 2 = 2, which in- 
dicates an excess of accretion in the equatorial region. We applied a mass 
threshold of 5 ■ IO^^A/q to our haloes for every redshift. Then, the halo 
population is different from one redshift to another. This selection effect 
may dominate the observed time evolution. 



3.3 One-point statistics: equatorial infall anisotropy 

The previous measurement doesn't account for dark matter falling 
into the halo with a very small angular momentum (radial orbits). 
We therefore measured the excess of equatorial accretion, Sm, de- 
fined as follow. We can measure the average flow density of matter, 
$r, in a ring centered on the equatorial plane: 



<E>r 



1 



pvr(r2) ■ dfl, 



(14) 



T/8<9-7r/2<7r/8 



where Sr 



dfl. The ring-region is defined by 



T/8<9-7r/2<7r/8 

the area where the polar angle satisfies 6poi ~ 7i"/2 ± n/8 which 
corresponds to about 40 % of the total covered solid angle. The 
larger this region is, the better the convergence of the average value 
of $r- , but the lower the effects of anisotropy, since averaging over 
a larger surface leads to a stronger smoothing of the field. This 
value of ±7r/8 is a compromise between these two contradictory 
trends. In the next section and in the appendix, we discuss more 
general filtering involving spherical harmonics which are related to 
the dynamical evolution of the inner component of the halo. We 
also measure the flow averaged on all the directions <E>: 



47r 



(15) 



Since we are interested in accretion, we computed $r and $ using 
only the infalling part of the density flux of matter, where pvr(f2) ■ 
dfl < 0, ignoring the outflows. We therefore define Sm as 



1 ■ 



This number quantifies the anisotropy of the infall. It is positive 
when infall is in excess in the galactic-equatorial plane, while for 
isotropic infall Sm = 0. The quantity Sm can be regarded as being 
the 'flux density' contrast of the infall of matter in the ring region 
(formally it is the centered top-hat-filtered mass flux density con- 
trast as shown in appendix ICH . This measurement, in contrast to 
those of the previous section does not rely on some knowledge on 
the inner region of the halo but only on the properties of the envi- 
ronment. 

Fig.|5|displays the normalized distribution of Sm measured for 
50 000 haloes with a mass in excess of 5 • 10^^ Mq and for different 
redshifts (z = 1.8, 1.5, 0.9, 0.3, 0.0). The possible values for Sm 
range between (5m —1 and Sm ^ 1-5. The average value {Sm) 
of the distributions is statistically larger than zero (see also Fig.|6j. 
Here ( ) stands for the statistical expectation, which in this pa- 
per is approximated by the arithmetic average over many haloes in 
our simulations. The antisymmetric part of the distribution of Sm is 
positive for positive Sm • The PDF of Sm is skewed, indicating an 
excess of accretion through the equatorial ring. The median value 
for Sm is (5mcd ~ 0.11, while the first 25% haloes have Sm < 
S25 = —0.11 and the first 75% haloes have Sm < S75 = 0.37. 
Therefore we have ((Jts — 5mod)/(525— <5mcd) = 1.13, which quan- 
tifies how the distribution of Sm is positively skewed. The skewness 
S3 = {{S - 5f)/{iS ~ is equal to 0.44. Combined with 

the fact that the average value {Sm) is always positive, this shows 
that the infall of matter is larger in the equatorial plane than in the 
other directions. 

This result is robust with respect to time evolution (see Fig. 
|6j. At redshift z — 1.8, we have {Sm) ~ 0.17 which falls down 
to {Sm) = 0.145 at redshift z — 0. This redshift evolution can be 
fitted as {Sm){z) = 0.0161(±0.0103)z + 0.147(±0.005). This 
trend should be taken with caution. For every redshift z we take 
in account haloes with a mass bigger than 5 ■ 10^^ Mq. Thus the 
population of haloes studied at z=0 is not exactly the same as the 
one studied at z=2. Actually, at z=0, there is a strong contribution 
of small haloes (i.e. with a mass close to 5 • IO^^A/q) which just 
crossed the mass threshold. The accretion on small haloes is more 
isotropic as shown in more details in appendix ID2I One possible 
explaination is that they experienced less interactions with their en- 
vironment and have had since time to relax which implies a smaller 
correlation with the spatial distribution of the infall. Also bigger 
haloes tend to lie in more coherent regions, corresponding to rare 
peaks, whereas smaller haloes are more evenly distributed. The 
measured time evolution of the anisotropy of the infall of matter 
therefore seems to result from a competition between the trend for 
haloes to become more symmetric and the bias corresponding to a 
fixed mass cut. 

In short, the infall of matter measured at the virial radius in 
the direction orthogonal to halo spin is larger than expected for an 
isotropic infall. 



3.4 Harmonic expansion of anisotropic infall 

As mentioned earlier (and demonstrated in appendix]^, the dy- 
namics of the inner halo and disk is partly governed by the statisti- 
cal properties of the flux densities at the boundary. Accounting for 
the gravitational perturbation and the infalling mass or momentum 
requires projecting the perturbation over a suitable basis such as the 
spherical harmonics: 



zu{fl) ^y^aTYrifl) 



(16) 



(17) 
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Figure 7. The convergence of modulus of the real part of (a™), for I = 
2,4,6, 8. The decomposition was computed for 25 000 haloes and each 
coefficient has been normalized with the corresponding Co (see text for 
details). Here, ( ) stands for the median while the error bars stand for the 
distance between the 5th and the 95th centile. The median value of (a™) is 
zero except for the (a^') coefficient: this is a signature of a field invariant to 
azimuthal rotations. 




7C/2 



37C/4 



7t 



Figure 8. An illustration of the convergence of the presented in Fig. 
171 The solid line stands for the azimuthal average of the spherical contrast 
of accretion computed using Eq. 1221 . The dotted line for the spheiical field 
reconstructed with £ ^ 5. The insert represents the reconstructed spherical 
field using the expansion of the ag^ of the mass flux measured at the virial 
sphere. The sphere presents an excess of accretion in the equatorial region 
because of the non-zero average value of chq coefficients (for even values 
of£). 



Here, vj stands for, e.g. the mass flux density, the advected mo- 
mentum flux density, or the potential perturbation. The resuhing 
a™ coefficients correspond to the spherical harmonic decomposi- 
tion in an arbitrary reference frame. The different m correspond 
to the different fundamental orientations for a given muhipole I. A 
spherical field with no particular orientation gives rise to a field av- 
eraged over the different realisations which appear as a monopole, 
i.e. (aim) = for ^ 7^ 0. Having constructed our virial sphere 
in a reference frame attached to the simulation box, we effectively 
performed a randomization of the spheres' orientation. However, 
since the direction of the halo's spin is associated to a general pre- 
ferred orientation for the infall, it should be traced through the a£m 
coefficients. Let us define the rotation matrix, TZ, which brings the 
z-axis of the simulation box along the direction of the halo's spin. 
The spherical harmonic d ecomposition centered on the spin of the 
halo, a™, is given by (e.g. lVarshalovich et"al]il988h ): 



7^: 



(18) 



If the direction of the spin defines a preferential plane of accretion, 
the corresponding aj^ will be systematically enhanced. We there- 
fore expect the equatorial direction (which corresponds to m = 
for every i) not to converge to zero. 

We computed the spherical harmonic decomposition of 
pVr{'&, (fi) given by Eq. <17t for the mass flux density of 25 000 
haloes at 2; = 0, up to £ = 15. For each spherical field of the mass 
density flux, we performed the rotation that brings the halo's spin 
along the z-direction to obtain a set of 'centered' aj^ coefficients. 
We also computed the related angular power spectrums C'e : 



Ce = 



4tv2£ + 1 ^ 



I = 



Ait 21 + 1 



E 



1 m 1 2 



(19) 



Let us define the normalized 
pendix IcTl . 



at 



(or harmonic contrast, see ap- 



(20) 



sign(ajjj 

This compensates for the variations induced by our range of masses 
for the halo. For each I, we present in Fig. the median value, 
\{Re{aT})\ for i = 2,4,6,8 computed for 25 000 haloes. Afl 
the 5™ have converged toward zero, except for the aio coefficients. 
The imaginary parts of have the same behaviour, except for 
the Im{5fo} coefficients which vanish by definition (not shown 
here). The ra = Q coefficients are statistically non-zero. We find 
{al) = -0.65±0.04, {ai) = 0.12±0.02, {&%) = -0.054±0.015 
and (ag) — 0.0145 ± 0.014. Errors stand for the distance between 
the 5th and the 95th centile. The typical pattern corresponding an 
m = Q harmonic is a series of rings parallel to the equatorial plane. 
This confirms that the accretion occurs preferentially in a plane per- 
pendicular to the direction halo's spin. 

The spherical accretion contrast {Spvr{^,<t>)) can be recon- 
structed using the (fi™) coefficients (as shown in the appendix): 



In Fig.|8| the polar profile: 

{5[,.^]W) = ^{ar)yr(^,o)-i, 



(21) 



(22) 



of this reconstructed spherical contrast is shown. This profile has 
been obtained using the (5™) coefficients with £ ^ 5 and £ ^ 15. 
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The contrast is large and positive near = 7r/2 as expected for an 
equatorial accretion. The profile reconstructed using ^ ^ 5 is quite 
similar to the one using I ^ 15. This indicates that most of the 
energy involved in the equatorial accretion is contained in a typical 
angular scale of 36 degrees (a scale which is significantly larger 
than 7r/20 corresponding to the cutoff frequency in our sampling 
of the sphere as mentionned earlier). 

Using a spherical harmonic expansion of the incoming mass 
flux density (Eq. <17H . we confirmed the excess of accretion in the 
equatorial plane found above. This similarity was expected since 
these two measurements (using a ring or using a spherical har- 
monic expansion) can be considered as two different filterings of 
the spherical accretion field as is demonstrated in the appendix Icl 
The main asset of the harmonic filtering resides in its relevance for 
the description of the inner dynamics as is discussed in section 6. 

3.5 Summary 

To sum up, the two measurements of section 3.2 and 3.3 (or 3.4) 
are not sensitive to the same effects. The first measurement (in- 
volving the angular momentum pL at the virial radius) is mostly a 
measure of the importance of infalling matter in building the halo's 
proper spin. The second and the third measurements (involving the 
excess of accretion in the equatorial plane, S,n, using rings and har- 
monic expansion) are quantitative measures of coplanar accretion. 
The equatorial plane of a halo is favoured relative to the accretion 
of matter (compared to an isotropic accretion) to a level of ~ 12% 
between z — 2 and z — 0. Down to the halo scale (~ 500 kpc), 
anisotropy is detected and is reflected in the spatial configuration 
of infalling matter. 



4 ANISOTROPIC INFALL OF SUBSTRUCTURES 

To confirm and assess the detected anisotropy of the matter in- 
fall on haloes in our simulations, let us now move on to a dis- 
crete framework and measure related quantities for satellites and 
substructures. In the hierarchical scenario, haloes are built up by 
successive mergers of smaller haloes. Thus if an anisotropy in the 
distribution of infalling matter is to be detected it seems reasonable 
that this anisotropy should also be detected in the distribution of 
satellites. The previous galactocentric approach for the mass flow 
does not discriminate between an infall of virialised objects and a 
diffuse material accretion and therefore is also sensitive to satellites 
merging: one would need to consider, say, the energy flux density. 
Ho wever, it is not c l ear if satellites are markers of the general infall 
and lVitvitska et all i2002h did not detect any anisotropy at a level 
greater than 20%. 

The detection of substructures and satellites is performed us- 
ing the code ADAPTAHOP which is described in details in the ap- 
pendix. This code outputs trees of substructures in our simulations, 
by analysing the properties of the local dark matter density in terms 
of peaks and saddle points. For each detected halo we can extract 
the whole hierarchy of subclumps or satellites and their characteris- 
tics. Here we consider the leaves of the trees, i.e. the most elemen- 
tary substructures the haloes contain. Each halo contains a 'core' 
which is the largest substructure in terms of particles number and 
'satellites' corresponding to the smaller ones. We call the ensemble 
core -I- satellites the 'mother' or the halo. Naturally the number of 
substructures is correlated with the mother's mass. The bigger the 
number of substructures, the bigger the total mass. Because the res- 
olution in mass of our simulations is limited, smaller haloes tend to 
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Figure 9. Excess probability, 1 + ^qs (Bcb), of the angle between the core's 
spin and the orbital momentum of satellites. Cores have at least 1 satel- 
lite (solid line), 4 satellites (dashed line) and 10 satellites (dotted line). 
These curves have been normalized by the expected isotropic distribution 
and the Gaussian correction was applied to account for errors on the angle 
determination. Here gcs {^) = is expected for an isotropic distribution 
of angles between the core 's spin and the orbital momentum. All satel- 
lites ai'e external to the core, yet an excess of alignment is present. The 
triangles represent the angle distribution, the error bars stand for the 3cr- 
dispersion for 50 subsamples of 10 000 satellites (out of 35 000) while the 
red curve stands for the best gaussian fit of the distribution for systems 
with at least one satellite (see Eg. 1 121 for parameterization). The best fit 
parameters are: ai = 0.3993 ± 0.0038, 02 = 0.0599 ± 0.0083, ag = 
0.8814 ± 0.0055, 04 = 0.9389 ± 0.0002. The isotropic case is excluded 
with a good confidence level, even for systems with a large number of satel- 
lites. 



have only one or two satellites. Thus in the following sections we 
will discriminate cases where the core have less than 4 satellites. 
A total of 50 000 haloes have been examined leading to a total of 
about 120 000 substructures. 



4.1 Core spin - satellite orbital momentum correlations 

In the mother-core-satellite picture, it is natural to regard the core 
as the central galactic system, while satellites are expected to join 
the halo from the intergalactic medium. One way to test the effect 
of large scale anisotropy is to directly compare the angle between 
the core's spin, Sc, and the satellites' angular momentum, Ls, rela- 
tive to the core. These two angular momenta are chosen since they 
should be less correlated with each other than e.g. the haloes' spin 
and the angular momentum of its substructures. Furthermore, par- 
ticles that belong to the cores are strictly distinct from those that 
belong to satellites, thus preventing any 'self contamination' ef- 
fect. As a final safeguard, we took in account only satellites with 
a distance relative to the core larger than the mother's radius. The 
core's spin is: 

Sc = ^(rp-rc) X (vp-Vc), (23) 

p 

where rp and Vp (resp. rc and Vc) stands for particles' position and 
velocities (resp. the core's centre of mass position and velocity) and 
where: 
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Tp < dc, 



(24) 



where dc is the core's radius. The angular momentum for a satelUte 
is computed likewise, with a different selection criterion on parti- 
cles, namely: 



(25) 



where stands for the satellite's centre of mass position and ds is 
its radius. 

Fig.|9|displays the reduced distribution of the angle, 6cb, be- 
tween the core's spin and the satellites' orbital momentum, where 
dcs is defined by: 



— 1/ Ls ■ Sc 



(26) 



The Gaussian correction was applied as described in section 3.2, to 
take into account the uncertainty on the determination of ^cs . 

The correlation of 6cs indicates a preference for the aligned 
configuration with an excess of ~ 12% of aligned configurations 
relative to the isotropic distribution. We ran Monte-Carlo realisa- 
tions using 50 subsamples of 10 000 haloes extracted from our 
whole set of substructures to constrain the error bars. We found 
a 3cr error of 6%: the detected anisotropy exceeds our errors, i.e. 
Ccs(^cs) is not uniform with a good confidence level. The varia- 
tions with the fragmentation level (i.e. the number of satellites per 
system) remains within the error bars. The best fit parameters for 
the measured distributions of systems with at least 1 satellite are 
ai = 0.3993 ± 0.0038, aa = 0.0599 ± 0.0083, ag = 0.8814 ± 
0.0055,04 = 0.9389 ± 0.0002 (see Eq.[T2lfor parametrization). 
Not surprisingly, a less structured system shows a stronger align- 
ment of its satellites' orbital momentum relative to the core's spin. 
In the extreme case of a binary system (one core plus a satellite), 
it is common for the two bodies to have similar mass. Since the 
two bodies are revolving around each other, a natural preferential 
plane appears. The core's spin will be likely to be orthogonal to 
this plane. Increasing the number of satellites increases the isotropy 
of the satellites' spatial distribution (the distributions maxima are 
lower and the slope toward low values of 9cs is gentler), but switch- 
ing from at least 4 satellites to at least 10 satellites per system does 
not change significantly the overwhole shape distribution. This sug- 
gests that convergence, relative to the number of satellites, has been 
reached for the 9cs distribution. 

As the measurements of the anisotropy factor 5m indirectly 
suggested, satellites have an anisotropic distribution of their direc- 
tions around haloes. Furthermore the previous analysis of the sta- 
tistical properties of 5 (section l33l indicated an excess of aligned 
configuration of 15% which is consistent with the current method 
using substructures. While the direction of the core's spin should 
not be influenced by the infall of matter, we still find the existence 
of a preferential plane for this infall, namely the core's equatorial 
plane. 



4.2 Satellite velocity - satellite spin correlation 

The previous sections compared haloes' properties with the prop- 
erties of satellites. In a galactocentric framework, the existence of 
this preferential plane could only be local. In the extreme each halo 
would then have its own preferential plane without any connection 
to the preferential plane of the next halo. Taking the satellite it- 
self as a reference, we have analyzed the correlation between the 
satellite average velocity in the core's rest frame and the struc- 
tures' spin. Since part of the properties of these two quantities are 



consequences of what happened outside the galactic system, the 
measurement of their alignment should provide information on the 
structuration on scales larger than the haloes scales, while sticking 
to a galactocentric point of view. 

For each satellite, we extract the angle, O^s, between the veloc- 
ity and the proper spin and derive its distribution using the Gaussian 
correction (see fig. llOt . The satellite's spin Ss is defined by: 



(rp-rs) X (vp-Vs 



(27) 



where Vs and Vs stands for the satellite's position and velocity in 
the halo core's rest frame. The angle, 9va, between the satellite's 
spin and the satellite's velocity is: 

Si, ■ Vj, , 



(28) 



Only satellites external to the mother's radius are considered 
while computing the distribution of angles. This leads to a sample 
of about 40 000 satellites, at redshift z ^ 0. The distribution ^{O^s) 
was calculated as sketched in section 2. An isotropic distribution of 
^vs would as usual lead to a uniform distribution £,{6vs) ~ 0. The 
result is shown in Fig. llOl The error bars were computed using the 
same Monte-Carlo simulations described before with 50 subsam- 
ples of 10 000 satellites. 

We obtain a peaked distribution with a maximum for O^s ~ 
■k/2 corresponding to an excess of orthogonal configuration of 5% 
compared to a random distribution of satellite spins relative to their 
velocities. The substructures' motion is preferentially perpendicu- 
lar to their spin. This distribution of angles for systems with at least 
1 satellite can be fitted by a Gaussian function with the follow- 
ing best fit parameters (see Ea ll2> : ai = 0.2953 ± 0.0040, ai — 
1.5447 ± 0.0015, aa = 0.8045 ± 0.0059, as = 0.9144 ± 0.0010. 
The variation with the mother's fragmentation level is within the 
error bars. However the effect of an accretion orthogonal to the di- 
rection of the spin is stronger for satellites which belong to less 
structured systems. This may be again related to the case where 
two comparable bodies revolve around each other, but from a satel- 
lite point-of-view. The satellite spin is likely to be orthogonal to the 
revolution plane and consequently to the velocity's direction. 

This result was already known for haloes in filaments 
(Faltenbacher et al. 1 2002)), where their motion occurs along the 
filaments with their spins pointing outwards. The current results 
show that the same behaviour is measured down to the satel- 
lite's scale. However this result should be taken with caution since 
Monte-Carlo tests suggest that the error (deduced from the 3a dis- 
persion) is about 4%. 

This configuration where the spins of haloes and satellites are 
orthogonal to their motion fit with the image of a flow of structures 
along the filaments. Larger structures are formed out of the merging 
of smaller ones in a hierarchical scenario. Such small substructures 
should have small relative velocities in order to eventually merge 
while spiraling towards each other The filaments correspond to re- 
gions where most of the flow is laminar, hence the merging between 
satellites is more likely to occur when one satellite catches up with 
another, while both satellites move along the filaments. During such 
encounter, shell crossin g induces vorticity perpendicu lar to the flow 
as was demonstrated in lPichon & BernardeatJ il999ft . This vortic- 
ity is then converted to momentum, with a spin orthogonal to the 
direction of the filament. 

Finally, the flow of matter along the filaments may also pro- 
vide an explanation for the excess of accretion through the equa- 
torial regions of the virial sphere. If a sphere is embedded in a 
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Figure 10. Excess probability. 1 + ^vs, of the angle between the substruc- 
tures' spin and their velocities in the mother's rest frame. The Gaussian cor- 
rection was applied to take into account uncertainty on the angle determina- 
tion. The distributions was measured for all mothers (solid), mothers with at 
least 4 substructures (dot) and mothers with at most 3 substructures (dash). 
The triangles represent the mean angle distribution. The error bars represent 
the Monte Carlo Scr-dispersion for 50 subsamples of 10 000 haloes (out of 
35 000). The red dash-dotted curve stands for the best fit of the distribution 
with a Gaussian function for systems with at least 1 satellite (see Ea. ll2l for 
parametrization). The best fit parameters are: a\ = 0.2953it0.0040, a\ = 
1.5447 ± 0.0015, 02 = 0.8045 ± 0.0059, 03 = 0.9144 ± 0.0010. In the 
core's rest frame, the satellites' motion is orthogonal to the direction of the 
satellites' spin. This configuration would fit in a picture where structures 
move along filamentary directions. 

'laminar' flow, the density flux detected near the poles should be 
smaller than that detected near the 'equator' of the sphere. The flux 
measured on the sphere is larger in regions where the normal to 
the surface is collinear with the 'laminar' flow, i.e. the 'equator'. 
On the other hand, a nil flux is expected near the poles since the 
vector normal to the surface is orthogonal to direction of the flow. 
The same effect is measured on Earth which receives the Sun radi- 
ance: the temperature is larger on the Tropics than near the poles. 
Our observed excess of accretion through the equatorial region sup- 
ports the idea of a filamentary flow orthogonal to the direction of 
the halo's spin down to scales ^ 500 kpc. 



5 PROJECTED ANISOTROPY 
5.1 Projected satellites population 

We looked directly into the spatial distributions of satellites sur- 
rounding the haloes cores to confirm the existence of a preferential 
plane for the satellites locations in projection. In Fig. ll II we show 
the compilation of the projected positions of satellites in the core's 
rest frame. The result is a synthetic galactic system with 100 000 
satellites in the same rest frame. We performed suitable rotations to 
bring the spin axis collinear to the z-axis for each system of satel- 
lites, then we added all these systems to obtain the actual synthetic 
halo with 100 000 satellites. The positions were normalized using 
the mother's radius (which is of the order of the virial radius). A 
quick analysis of the isocontours of the satellite distributions in- 
dicates that satellites are more likely to be found in the equatorial 



plane, even in projection. The axis ratio measured at one mother's 
radius is e{Rm) = a/b — 1 = 0.1 with a > h. We compared 
this distribution to an isotropic 'reference' distribution of satellites 
surrounding the core. This reference distribution has the same av- 
erage radial profile as the measured satellite distributions but with 
isotropically distributed directions. The result of the substraction of 
the two profiles is also shown in Fig. ll II The equatorial plane (per- 
pendicular to the z-axis) presents an excess in the number of satel- 
lites (light regions). Meanwhile, there is a lack of satellites along 
the spin direction (dark regions). This confirms our earlier results 
obtained using the alignment of orbital momentum of satellites with 
the core's spin, i.e. satellites lie more likely in the plane orthogonal 
to the halo spin d irection. Qualitatively, these results have already 
been obtained bv iTormei] <1997.) . where the major axis of the el- 
lipsoid defined by the satellite's distribution is found to be aligned 
with the cluster's major axis. This synthetic halo is more directly 
comparable to observables since, unlike the dark matter halo itself, 
the satellites should emit light. Even though ACDM predicts too 
many satellites, its relative geometrical distribution might still be 
correct. In the following sections, our intent is to quantify more 
precisely this effect. 

The propension of satellites to lie in the plane orthogonal to 
the direction of th e core's spin appears as an 'anti-Holmberg' effect. 
iHolmberd {l9l4\ and more recently IZaritskv et al.l i 19971) have 
found observationnally that the distribution of satellites around 
disks is biased towards the pole regions. Thus if the orbital mo- 
mentum vector of galaxies is aligned with the spin of their parent 
haloes, our result seems to contradict these observations. One may 
argue that satellites are easier to detect out of the galactic plane. 
Furthermore our measurements are carried far from the dis k while 
its influence is not taken in account. fHuang & Carlber j |'l997) have 
shown that the orbital decay and the disruption of satellites are 
more efficient for coplanar orbits near the disk. It would explain the 
lack of satellites in the disk plane. Thus our distribution of satellites 
can still be made consistent with the 'Holmberg effect'. 

5.2 Projected satellite orientation and spin 

In addition to the known alignment on large scales, we have shown 
that the orientation of structures on smaller scales should be dif- 
ferent from the one expected for a random distribution of orienta- 
tions. Can this phenomenon be observed ? The previous measure- 
ments were carried in 3D while this latter type of observations is 
performed in projection on the sky. The projection 'dilutes' the 
anisotropy effects detected using three-dimensional information. 
Thus an effect of 15% may be lowered to a few percents by project- 
ing on the sky. However, even if the deviation from isotropy is as 
important as a few percents, as we will suggest, this should be rel- 
evant for measurements involved in extracting a signal just above 
the noise level, such as weak lensing. 

To see the effect of projection on our previous measurements, 
we proceed in two steps. First, every mother (halo core -1- satellites) 
is rotated to bring the direction of the core's spin to the z-axis. Sec- 
ond, every quantity is computed using only the y and z components 
of the relevant vectors, corresponding to a projection along the x- 
axis. 

The first projected measurement involves the orientation of 
satellites relative to their position in the core's rest frame. The spin 
of a halo is statistically orthogo nal to the main axis of the distri- 
bution of matter of that halo I Faltenbacher et alj i2002h '). and as- 
suming that this property is preserved for satellites, their spin Sg 
is an indicator of their orientation. The angle, 6p (in projection). 
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Figure 11. The projected distribution of satellites around the core's centre 
of mass. We used the position of 40 000 satellites around their respective 
core to produce a synthetic halo+satellites (a 'mother') system. The pro- 
jection is performed along the x-axis. The y and z coordinates are given in 
units of the mother's radius. The z-axis is colinear to the direction of the 
core's spin. Top: The isocontours of the number density of satellites around 
the core's centre of mass present a flattened shape. The number of satel- 
lites is lower in darker bins than in Hghter bins. The flattened isocontours 
indicate that satellites he preferentially in the plane orthogonal to direction 
of the spin. Bottom: The excess number of satellites surrounding the core. 
We compared the distribution of satellites measured in our simulations to 
an isotropic distribution of satellites. Light zones stand for an excess of 
satellites in these regions (compared to an isotropic distiibution) while dark 
zones stand for a lack of satellites. The satellites are more numerous in 
the equatorial region than expected in an isotropic distribution of satellites 
around the core. Also, there are fewer satellites along the spin's axis than 
expected for an isotropic distribution of satellites. 



between the satellites' spin and their position vector (in the core's 
rest frame) is computed as follow: 



5p = cos (t^^ 



with 



(29) 



Vac = ra — Tc, (30) 

where and Tc stand respectively for the position vector of the 
satellite and the core's centre of mass. Two extreme situations 
can be imagined. The 'radial' configuration corresponds to a case 
where the satellite's main axis is aligned with the radius joining the 
core's centre of mass to the satellite centreof mass (spin perpendic- 
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Figure 12. Excess probability, 1 + gp, of the projected angles between the 
direction of the spin of substructures and their position vector in the core's 
rest frame . The projection is made along the x-axis where the z-axis is con- 
current with the core's spin direction. The solid line represents the average 
distribution of projected angles of 50 subsamples of 50 000 substructures 
(out of 100 000 available substructures). The error bars represents the 3cr 
-dispersion relative to these 50 subsamples. An isotropic distribution of ori- 
entation would con'espond to a value of 1 for 1 + ^p. The projection plus 
the reference to the position vector instead of the velocity's direction low- 
ers the anisotropy effect. The red curve stands for the best Gaussian fit of 
the excess probability (see Eg. 1 121 for parameterization). The best fit pa- 
rameters are: ai = 0.0999 ± 0.0030, a2 = 1.5488 ± 0.0031, a-j = 
0.8259 ± 0.0131, a4 = 0.9737 ± 0.0007. It seems that on average the 
projected orientation of a substructure is orthogonal to its projected posi- 
tion vector. 



ular to the radius, or 9p ~ 7r/2). The 'circular' configuration is the 
case where the satellite main axis is orthogonal to the radius (spin 
parallel to the radius, ^p ~ 0[7r]). These reference configuration 
will be discussed in what follows. 

The resulting distribution, 1 + ^p(6'p), is shown in Fie. 1121 
As before, an isotropic distribution of orientations would lead to 
^p(Sp) — 0. The distribution is computed with 100 000 satellites, 
without the cores, while the error bars result from Monte-Carlo 
simulations on 50 subsamples of 50 000 satellites each. As com- 
pared to the distribution expected for random orientations, the or- 
thogonal configuration is present in excess of ^p(7r/2) ~ 0.02. 
If the spin of satellites is orthogonal to their principal axis, the 
direction vector in the core's rest frame is more aligned with the 
satellites principal axes than one would expect for an isotropic 
distribution of satellites' orientations. This configuration is 'ra- 
dial'. The peak of the distribution is slightly above the error bars: 
A^p(6'p ~ 7r/2) ~ 0.02. The distribution can be fitted by the 
Gaussian function given in Eq. 1121 with the following parame- 
ters: ai = 0.0999 ± 0.0030, 02 = 1.5488 ± 0.0031, aa = 
0.8259 ± 0.0131, 04 = 0.9737 ± 0.0007. The alignment seems 
to be difficult to detect in projection. With 50 000 satellites, we 
barely detect the enhancement of the orthogonal configuration at 
the 3-(7 level, thus we do not expect a detection of this effect at the 
l-(7 level for less than 6000 satellites. Nevertheless, the distribution 
of the satellites' orientation in projection seems to be 'radial' on 
dynamical grounds, without reference to a tensing potential. 

Our previous measurement was 'global' since it does not take 
into account the possible change of orientation with the relative 
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Figure 13. Radial and azimuthal grid of the excess probability, 1 + ^p, 
of the projected angles between the direction of the spin of substructures 
and their direction relative to the central position of the core (as shown on 
average in Fig |12l . The projection is made along the x-axis where the z- 
axis is concurrent with the direction of the core's spin. Each row represents 
a distance relative to central core in the mother's radius units (from bottom 
to top): R e [0,0.4[, _R 6 [0.4,0.8[, R S [0.8, 1.2[, R e [1.2, 1.6[ 
and R S [1.6, 2[. Each column represents an angular distance (in degrees) 
relative to the direction of the core's spin (z-axis) : <f>s S]0, 36], <j>s G 
[mj2[. (t>s e [72, 108[, </is e [lOS, 144[ and S [144, 180[. The 
isotropic orientation distribution corresponds to a value of 1. Each sector 
presents a preferential direction that depends on its position relative to the 
spin direction of the central core. The distributions are computed using 50 
samples of 50 000 satellites each. In each sector, the points represents the 
distribution averaged over the 50 samples. The eiTor bars represent the 3(t- 
Monte Carlo dispersion of the distribution over these 50 samples. 

position of the satellites in the core's rest frame. In Fig. 1131 we 
explore the evolution of 1 + with the radial distance relative 
to the core's centre of mass and with the angular distance relative 
to the z-axis, i.e. relative to the direction of the core's spin. The 
previous synthetic halo was divided in sectors and for each sector, 
1+^p can be computed. The sectors are thus defined by their radius 
(in the mother's radius units): R ^ 0.4, 0.4 < J? < 0.8, 0.8 < 
R ^ 1.2, 1.2 < 7? ^ 1.6 and 1.6 < < 2 and by their polar 
angle relative to the direction of the core's spin (in degrees): (f)s ^ 
36, 36 < 0^ < 72, 72 < (jjs ^ 108, 108 < (t>s 144 and 
144 < (j)a ^ 180. Each of the previous Monte-Carlo subsamples 
can also be divided into sectors in order to compute the dispersion 
cr for the distributions within the subsamples. The error bars still 
represent the 3(j dispersions. 

The Fie. ll4l is a qualitative representation of the results pre- 
sented in Fig. 1131 Each sector with _R ^ 1 in Fig. ll3l is represented 
by an ellipse at its actual position. The orientation of the ellipse is 
given by the angle of the maximum of the corresponding 1 +5p {Op) 
function. We chose to represent the spin's direction perpendicular 
to the ellipse's major axis. We also chose to scale the ellipse axis 
ratio with the signal-to-noise ratio of 1 + ^p(Sp). Indeed large er- 
rors leads to weak constrains on the spin orientation and the galaxy 
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Figure 14. Geometric configuration of mean satellites around their core 
galaxy; each panel of Fig. ll3l is represented by an ellipse at its log radius and 
angle around the core galaxy. The axis ratio of the ellipses is proportional 
to the peak-to-peak amplitude of the corresponding correlation (accounting 
for the relative SNR), while its orientation is given by the orientation of the 
maximum of 1 + 5p- 



would be seen as circular on average. Conversely a strongly con- 
strained orientation leads to a typical axis ratio of 0.5. 

Two effects seem to emerge from this investigation. For some 
sectors, the orthogonal configuration is in excess compared to an 
isotropic distribution of satellites' orientation relative to the radial 
vector. This seems to be true especially for radii smaller than the 
mother's radius but the effect is still present at larger distances, es- 
pecially near <j)s ~ 7r/2. Switching from low values to high values 
of (j)a changes the slope of the 1 +^p {6) distribution. This may be a 
marker of a 'circular configuration' of the orientation of satellites. 

The existence of a 'radial' component in the orientation of the 
satellites was expected, both from the unprojected measurements 
made in the previous sections and from the global distribution ex- 
tracted from the projected data. The fact that the 'radial' signature 
is stronger around the equatorial plane (72 < (f>s < 108 in Fig. 
I13> may be an another evidence for a filamentary flow of satellites, 
even in projection. It seems that the existence of a 'circular' compo- 
nent was mostly hidden in the previous measurements by the dom- 
inant signature of the 'radial' flow. Nevertheless, the dominance of 
'circular' orientations near the poles fits with the picture of a halo 
surrounded by satellites with their spin pointing orthogonally to the 
filament directions. 

The 'circular' flow may alternatively be related to the flow of 
structures around clusters located at the connection between fila- 
ments. There are observations of such configurations (Kitzbichler 
& Saurer 2003), where galaxies have their spin pointing along their 
direction of accretion and these observations could be consistent 
with our 'circular' component. 



6 APPLICATIONS 

Let us give here a quick overview of the implications of the pre- 
vious measurements for the inner dynamics of the halo down to 
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galactic scales. In particular let us see how the self-consistent dy- 
namical response of the halo propagates anisotropic infall inwards, 
and then briefly and qualitatively discuss implications of anisotropy 
to galactic warps, disk thickening and lensing. 



6.1 Linear response of galaxies 

In the spirit of e.g. lKalnaisI h97lh or lTremaine & WeinberEl<I984 
we show in appendix]^ and elsewhere ( Aube rt & PichonI yOOJ)) 
how to propagate dynamically the perturbation from the virial ra- 
dius into the core of the galaxy using a self consistent combination 
of the linearized Boltzmann and Poisson equations under the as- 
sumption that the mass of the perturbation is small compared to the 
mass of the host galaxy. Formally, we have: 



r(x, t) = R[F, n,x,t- t] r)) 



(31) 



where R is a linear operator which depends on the equilibrium state 
of the galactic halo (-l-disk) characterized by its distribution func- 
tion F, and r(x, t) represents the self consistent response of the 
inner halo at time t due to a perturbation zu{fl, r) occuring at time 
r. Here -ou represents formally the perturbed potential on the virial 
sphere and the flux density of advected momentum, mass and ki- 
netic energy at i?200- A 'simple' expression for R is given in Ap- 
pendix|^for the self consistent polarisation of the halo. The linear 
operator, R, follows from Eq. <A6> . <A13> and <A16t . These equa- 
tions generalize the work of Kalnajs in that it accounts for a con- 
sistent infall of advected quantities at the outer edge of the halo. It 
is shown in particular in appendix]^ that self-consistency requires 
the knowledge of all ten (scalar, vector and symmetric tensor) fields 
TOp(fi, r), 7J7pv(n, r), tUpCT^CT^ (n, r). 

When dealing with disk broadening, r could be the velocity 
orthogonal to the plane of the disk, or, for the warp, its amplitude, 
as a function of the position in the disk, x (or the orientation of 
each ring if the warp is described as concentric rings). More gen- 
erally, it could correspond to the perturbed distribution function of 
the disk-l-halo. The whole statistics of r is relevant. The average 
response {r(x, t)) can be written as: 



{r(x, t)) = R(t^(n, r)) = ^YrmiaT) 



(32) 



Since the accretion is anisotropic, (a™) do not converge toward 
zero (see section l7!4l inducing a non-zero average response. Most 
importantly the two point correlation of the response since it will 
tell us qualitatively what the correlation length and the root mean 
square amplitude of the response will be. For the purpose of this 
section, and to keep things simple, we will ignore temporal issues 
(discussed in Appendix IaI altogether, both for the mean field and 
the cross-correlations. The two-point correlation of r(x) then de- 
pends linearly on the two-point correlation of zo: 



{r(x) ■ r^(y)) = R(i:^(fl) ■ 



(33) 



where stands for the transposition. Clearly, if the infall, ii7(f2), 
is anisotropic the response will be anisotropic. As was discussed in 
section 3.4 when the infall is not isotropic, we have 
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Figure 15. The residual anisotropic harmonic power spectra, A.R^", intro- 
duced in Eq. 1351 as a function of m for (. = 1, 2, 3, 4. These residuals will 
serve as input to the computation of the dynamical response of the halo. 



which would be identically zero if the field were stationary on the 
sphere. Here AT?™ represents the anisotropic excess for each har- 
monic correlation. In particular, the excess polarisation of the re- 
sponse induced by the anisotropy reads 



A{r(x).r^(y)) 



(36) 



Fig. fTsl displays ARf, for I = 1,2,3,4. The different AEf" 
clearly converge toward different non-zero values. Consequently 
the response should reflect the anisotropic nature of the external 
perturbations. 

It is beyond the scope of this paper to pursue the quantitative 
exploration of the response of the inner halo to a given anisotropic 
infall, since this would require an explicit expression of the re- 
sponse operator, R, for each dynamical problem investigated. 

6.2 Implication for warps, thick disks and lensing 

In this paper, the main emphasis is on measured anisotropics. It 
turns out that it never exceeds 15 % in accretion. For a whole class 
of dynamical problems where anisotropy is not the dominant driv- 
ing force it can be ignored at that level. Here we now discuss qual- 
itatively the implication of the previous measurements to galactic 
warps, the thick disk and weak lensing where anisotropy is essen- 
tial. 



6.2.1 Galactic warps 

The action of the torque applied on the disk of a galaxy is dif- 
ferent for different angular and radial position of the perturbation. 
Consequently the warp's orientation and its amplitude are functions 
of the spatial configuration of the external potential. For example, 
iLopez-Corredoira et alji20o3) found that the warp's amplitude due 
to an intergalactic flow is dependent on the direction of the incom- 
ing 'beam' of matter. Having modelled the intergalactic flow ap- 
plied to the Milky Way, they found that the warp amplitude rises 
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steeply as the beam leaves the region coplanar to the disk and this 
warp amplitude reaches a maximum for an inclination of 30 de- 
grees relative to the disk's plane. As the beam direction becomes 
perpendicular to the galactic plane, the warp amplitude decreases 
slowly. In this context, the existence of a typical spatial configu- 
ration for the incoming intergalactic matter or satellites infall may 
induce a kind of 'typical' warp in the disk of galaxies. 

The existence of a preferential plane for the accretion of angu- 
lar momentum also implies that the recent evolu tion of the halo's 
spin has been rather smooth, feullock et alj i200ll) have shown that 
the angular momentum tends to remain aligned within haloes. Fur- 
thermore, the accretion of matter by haloes is preferentially per- 
formed on plunging radial orbits, thus the inner parts of haloes are 
aware of the properties of the recently accreted angular momen- 
tum. Therefore, a disk embedded in the h a lo wo uld also 'feel' this 
anisotropic accretion. lOstriker & Binned il989 h have shown that 
the misalignment of the accreted angular momentum and the disk's 
spin forces the latter to slew the symmetry axis of its inner parts. 
The warp line of nodes is also found to be align ed with t he axi s 
of the torque applied to the disk. As stressed by iBinnevI il992h . 
non straight line of nodes can be associated with changes in the 
direction of the acc reted angular momentum. Using a sample of 
12 galaxies. iBriggj ( 1990) established rules of thumb for galactic 
warps, one of them being that the line of nodes is straight in the in- 
ner region of disk while it is wound in the outer parts. If the angular 
momentum is accreted along a stationary preferential direction, as 
we suggest, the warp line of nodes should remain mostly straight. 
However, if the accretion plane differs slightly from the disk plane, 
more than one direction of accretion become possible (by sym- 
metry around the vector defining the disk plane) and, as a conse- 
quence, different directions are possible for the torque induced by 
accreted matter. We may then consider a varying torque along ac- 
cretion history with an accreted angular momentum 'precessing' 
around the halo's spin but close to its direction. In this scenario, the 
difference in the behaviour of the warp line of nodes between the 
inner and outer regions of the galaxies may be explained. 

6.2.2 Galactic disk thickening 

Thin galactic disks put serious constraints on merging scenarii, 
since their presence implies a fine-tuning between the cooling 
mechanisms {e.g. coplanar infall of gas), and the heating processes 
(merging of small virialised objects, deflection of spirals on molec- 
ular clouds). It has been shown that sma ll merger s can produce a 
thick disk (e.g. .Ouinn et al. (.1993.,) , Walker et alJ ri996) ). How- 
ever, the presence of old stars within the thin disk cannot be ex- 
plained in the framework of the merging scenario unless a fraction 
of the accretion took place within the equatorial plane of the galaxy. 
Furthermore, the geometric characteris tic of the infall is essential i n 
the formation process of a thick disk. In lVelazguez & White! il999l) . 
numerical simulations of interactions between galactic disks and 
infalling satellites show that the heating and thickening is more ef- 
ficient for coplanar satellites. They also stressed the differences be- 
tween the effect of prograde or retrograde orbits of infalling satel- 
lites (relative to the rotation of the disk): prograde orbits induce disk 
heating while retrograde orbits induce disk tilting. Our results indi- 
cate that the infall is preferentially prograde and coplanar relative to 
the halo's spin: if we consider an alignment between the halo's spin 
and the galaxy's angular momentum, the thickening process may be 
more efficient than the one expected in an isotropic configuration of 
infalling matter. Furthermore, our estimate of the fraction of copla- 
nar accretion at the virial scale may be considered as a lower bound 



near the disk since the presence of a disk will focus the infall closer 
to the galactic plane. In fact, Huang & Carlberg ( 1997) found that 
the disk tends to tilt toward the orbital plane of infalling prograde 
low-density satellites. This effect would also contribute to enhance 
the excess of coplanar accretion down to galactic scales. 

However the nature of infalling virialised objects was shown 
to affe ct their ability to heat or destroy the disk. Huang & Carlbe^ 
( 1997) found that the presence of low density satellites should in- 
duce preferentially a tilting of the disk instead of a thickening: one 
needs to enhance the relative mass of the satellite (~ 30% of the 
disk mass) to produce an observable thickening in the inner parts 
of the galaxy. Unfortunately such a massive satellite has a destruc- 
tive impact on the outer parts of the disk. The relationship between 
the excess of accretion and the satellite mass should be constrained 
but our limited mass resolution prevents us from performing such a 
quantitative analysis. We should therefore aim at achieving higher 
angular resolution of the virial sphere and higher mass resolution 
in order to describe well compact virialised objects. 



6.2.3 Gravitational lensing 

The first detection of cosmic shea r was r eported by f o ur dif - 
ferent groups in 2000 ('Bacon et al.' !200(f), 'Kaiser et al J llOOtj) . 
Van Waerbeke et al. ( 2OO0), Wittman et al. ( 2000)). One of the ba- 
sic assumptions made by cosmic shear studies is that the intrin- 
sic ellipticities of galaxies are expected to be uncorrelated, and 
that the observed correlations are the results of gravitational lens- 
ing induced by the large scale structures between those galaxies 
and the observer. Hence, the detection of weak lensing signal as- 
sumes a gravitationally induced departure from a random distri- 
bution of the galactic shapes. Consequently, if there exists intrin- 
sic alignments or preferential patterns in galactic orientations, this 
would potentially affect the interpretation from weak lensing mea- 
surements. Several papers have already considered the 'contamina- 
tion' of the weak lensin g signal by int r insic g alactic alignment. Us- 
ing analytic arguments. ICatelan et alj i200ll) have shown that such 
alignments should exist. The issue of the amplitude of the intrin- 
sic correlations compared to the correlation induced by the cos- 
mic shear has also b een explored by .Croft & M etzlet (2000) and 
iHeavens et alj J2000l) . The 'intrinsic' correlations may overcome 
the shear-induced signal in surveys with a narrow redshift range. 
We have shown that the orientation of satellites around haloes is 
not randomly distributed, which is a clear indication of intrinsic 
correlations for our considered scales (~ 500 kpc). Taking z,n = 1 
as a typical median redshift for large lensing surveys, the corre- 
sponding angular scale is 1 arcminute in our simulations' cos- 
mogony. Furthermore, the prospect of studying the redshift evo- 
lution of gravitational clustering via shear measurements will re- 
quire investigating narrower redshift bins and as such, small scale 
dynamically induced polarisation m ight become an issue. As rec- 
ommended by Catela n et alj i200 lh. our measurement may also be 
used as a 'numerical' calibration of the relation between ellipticity 
and tidal fields. Interestingly, they suggested to compensate for the 
finite number of galaxies around clusters by 'stacking' several clus- 
ters, which is precisely the procedure we followed to extract signal 
from our simulations. Finally, Weak lensing predicts no 'curl' com- 
ponent in the shear field (e.g. Pen et al. ( 2000)) and such 'curl' con- 
figurations would serve to extract the intrinsic signal. Even though 
satellites exhibit both 'circular' and 'radial' configurations in our 
simulations, we do not observe a clear signature of a 'curl' compo- 
nent of orientations at our level of detection. 



Dark matter anisotropic cosmic infall on haloes 15 



Virial orbital 
momentum 




Core spin & orbital momentum 
at virial radius are aligned 



Satellite 



orbital \ j 
momentum II Core Spin 



oogo^oi 



o 



o 



o 



Core Spin & satellites orbital 
momentum aligned 




Excess of accretion in equatorial 
plane (ring & harmonic ) 




Satellites motion is perpendicular 
to spin direction 



1 1 Core Spin 




In projection, satellites 
lies in plane orthogonal 
to core spin 



^ Core Spin^ 



cPo" 



o 

'Radial' orientation in 
equatorial plane. 
'Circular' orientation near 
poles 



Figure 16. A schematic representation of all estimate of anisotropic accretion considered in this paper. 1): we measured the distribution of the angle between 
the orbital momentum on the virial sphere and the halo's spin. The average orbital momentum measured on the virial sphere is mostly aligned with the spin of 
the halo embedded in the virial sphere (discussed in section lTll . 2): we compared the accretion in the plane orthogonal to the direction of the halo's spin with 
the average accretion on the sphere. On the virial sphere, we detected an excess of ring-like or harmonic accretion in the equatorial plane (discussed in section 
I3.3l and l3!4l . 3): in projection, we used a 'synthetic' halo to look at the distribution of satellites detected with ADAPTAHOP and at the orientation of their spin 
around the direction of the spin. In projection, satellites lie preferentially in the projected equatorial plane (discussed in section lsTI . 4): We measured the angle 
between the halo's spin and the orbital momentum of each satellite. The orbital momentum of satellites is preferentially aligned with the spin of their hosting 
core (discussed in section l4m . 5): We compared the orientation of each satellite velocity vector (in the core's rest frame) with the orientation of their own 
spin. The velocity vector of satellites (in the core's rest frame) is orthogonal to the direction of their spin (discussed in section l4!2l . 6): In the equatorial plane, 
the projected orientation of satellites is more 'radial', while near the direction of the spin a 'circular' configuration of orientation seems to emerge (discussed 
in section ls^ . 



Table 1. Summary of the fitting parameters for the angular correlations 
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2.351 ±0.006 


-0.178±0.002 
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0.728±0.001 
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0.399±0.003 


0.059±0.008 


0.881 ±0.005 


0.938±0.000 


Ovs 


0.295±0.004 


1.544±0.001 


0.804±0.005 


0.914±0.001 


Op 


0.099 ± 0.003 


1.548 ±0.003 


0.825 ± 0.013 


0.973 ± 0.000 



Here SpL is the angle between the halo's spin and the angular momentum 
measured on the virial sphere; 9pvrL is the angle between the halo's spin 
and the accreted angular momentum measured on the virial sphere; Bcs is 
the angle between the core's spin and the satellite orbital momentum; Ovs 
is the angle between the satellite velocity in the core 's rest frame and the 
satellite's spin; Op is the projected angle between the satellite's spin and 
its direction relative to the core s' position. The fitting model we used is 
1 + m = ai/(v^a3) exp [-(0 - a2)V{2a^)] + 14- 



Table 2. Summary of other quantities related to anisotropic accretion 





0.0161(± 0.0103) z +0.147 (± 0.005) 




0.44 


e{Rm) 


0.1 


a20 


-0.65 ± 0.04 


5,4,0 


0.12 ± 0.02 


0,60 


-0.054 ±0.015 


0,80 


0.0145 ± 0.0014 



(Sm) (z) is the redshift evolution of the average excess of accretion in the 
plane orthogonal to the direction of the spin. S3 ((5m ) is the skewness of the 
distribution of excess of accretion. e{Rm) is the axis ratio a/b — 1 with 
a > b of the projected satellite distribution. 020 , <i40 > 060 and aso are the 
normalized harmonic coefficients of the 'equatorial' modes. 

7 CONCLUSION & PROSPECTS 



7.1 Conclusion 



Using a set of 500 ACDM simulations, we investigated the proper- 
ties of the spatial configuration of the cosmic infall of dark matter 
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around galactic ~ L* haloes. The aim of the present work was to 
find out if the existence of preferential directions existing on large 
scales (such as filaments) is reflected in the behaviour of matter 
accreted by haloes, and the answer is a clear quantitative yes. 

Two important assumptions were made in the present pa- 
per. We did not consider different class of haloes' masses (except for 
Fie. lD2> . but instead applied normalisations to includes all haloes in 
our measurements (considering e.g. the statistical average of con- 
strasts). We also did not take in account outflows and focused on 
accreted quantities. 

First we looked at the angular distribution of matter at the in- 
terface between the intergalactic medium and the inner regions of 
the haloes. We measured the accreted mass and the accreted an- 
gular momentum at the virial radius, describing these quantities as 
spherical fields. 

• The total (resp. advected) angular momentum measured at the 
virial radius is strongly aligned with the inner spin of the halo with 
a proportion of aligned configuration 30% (resp. 50%) more fre- 
quent than the one expected in an isotropic distribution of accreted 
angular momentum (1 + 5ls(0) ~ 1.5). This result reflects the 
importance of accreted angular momentum in the building of the 
haloes' inner spin. 

• The accretion of mass measured at the virial radius in the ring- 
like region perpendicular to the direction of the halo's spin is ~ 
15% larger than the one expected in the case of an isotropic infall 
of matter. 

We also detected the excess of accretion at the same level in the 
equatorial plane using a spherical harmonic expansion of the mass 
density flux. 

• In the spin's frame, the average of the harmonic a^o coeffi- 
cients does not converge toward zero, indicating that there is a sys- 
tematic accretion structured in rings parallel to the equatorial plane. 

Using the substructure detection code ADAPTAHOP, we con- 
firmed that the existence of a preferential plane for the infalling 
mass is reflected in the distribution of satellites around haloes. 

• Investigating the degree of aligrmient between the orbital mo- 
mentum of satellites and the central spin of the halo, it is shown that 
the aligned configuration is present in excess of ~ 12%. Satellites 
tend to revolve in the plane orthogonal to the direction of the halo's 
spin. The two methods (using spherical fields and satellites' detec- 
tion) yield consistent results and suggest that the image of a spher- 
ical infall on haloes should be reconsidered at the quoted level. 
We studied the distribution of the angle between the direction of 
accretion of satellites and their own spin. 

• An orthogonal configuration is 5% more frequent than what 
would one expect for an isotropic distribution of spin and direc- 
tions of accretion. Satellites tend to be accreted in the direction 
orthogonal to their own spin. 

These findings are interpreted as the results of the filamentary 
flows of structures, where satellites and haloes are accreted along 
the main direction of filaments with their spins orthogonal to this 
preferential direction. The flow along filaments also explains why 
the matter is accreted preferentially in the equatorial plane at the 
virial radius. The halo points its spin perpendicular to the flow and 
sees a larger flux in the regions normal to the flow direction, i.e. 
near the equator. Thus, it appears that the existence of preferen- 
tial directions on large scales is still relevant on galactic scales and 
should have consequences for the inner dynamics of the halo. 

We addressed the issue of observing these alignments in projec- 
tion. 



• The distribution of satellites projected onto the sky is flattened, 
with an axis ratio of 1.1 at the virial radius. 

• It seems that the orientation of satellites around their haloes 
is not random, even if the two dimensional representation dilutes 
the effects of alignments. The 'radial' orientation, where the satel- 
lites main axis is aligned with the line joining the satellite to the 
halo centre, is ~ 5% more frequent than the one expected in a 
completely random distribution of orientation. The 'circular' con- 
figuration, where the satellites main axis is perpendicular to the line 
joining the satellite to the halo centre, is also present in excess com- 
pared to an random distribution near the pole of the host galaxy. 

All corresponding fits are summarized in Table 1 and 2, while Fig. 
1 161 gives a schematic view of the measurements we carried out. 

We investigated how the self-consistent dynamical response 
of the halo would propagate anisotropic infall down to galac- 
tic scales. In particular we gave the corresponding polarisation 
operator in the context of an opened system. We have shown 
in appendix |A| that accounting for dark matter infall required 
the knowledge of the first three moments of the flux densities, 

■CUp{Q,, r), tI7pv(0, T),U]p„^a. (fl, r). 

It is suggested that the existence of a preferential plane of ac- 
cretion of matter, and thus of angular momentum, should have an 
influence on warp generation and disk thickening. If the anisotropic 
properties of infalling matter measured in the outer parts of haloes 
are conserved in the inner region of galaxies, there may exist a 'typ- 
ical' warp amplitude and this anisotropic accretion of matter may 
explain the properties of warp line of nodes. In the same spirit, the 
efficiency of the thickening of the disk may be enhanced or reduced 
by equatorial accretion. Finally, our finding of intrinsic alignments 
on small scales as well as specific orientations of structures should 
be relevant for cosmic shear studies on wide and shallow surveys. 

7.2 Prospects 

The main purpose of our investigation was to provide quantitative 
measurements of the level of anisotropy involved in the infall on 
scales 500 kpc. The next step should clearly involve working 
out quantitatively their implications for warp, disk heating etc... as 
was discussed in section|6| 

Our measurements were carried out at R200, which on galac- 
tic scales is a long way from the inner region of the galaxy. One 
should clearly propagate the infall (and its anisotropy) towards the 
centre of the galaxy, and more radial infalling components will 
play a more important role and should be weighted accordingly. 
It should also be stressed that we did not take into account the extra 
polarisation induced by the presence of an embedded disk, which 
will undoubtedly reinforce the polarisation and the anisotropy of 
the infall. We also concentrated on mass accretion, as the lowest 
order moment of the underlying "fluid" dynamics. Clearly higher 
moments involving the anisotropically accreted momentum, the ki- 
netic energy etc. are dynamically relevant for the evolution of the 
central object as is discussed in section 6 and in the appendix. The 
time evolution of the statistics of these flux densities is also es- 
sential for the inner dynamic of the halo and should be addressed 
systematically as well. It will be worthwhile to explore different 
cosmologies and their implications on small scale dynamics, and 
on the characteristics of infalling clumps, though we hope that the 
qualitative results sketched here should persist. 

It should be emphasized that some aspects of the present work 
are exploratory only, in that the resolution achieved (Mhaio > 
5 ■ 10^^ Mq) is somewhat high for galaxies. In fact, it would 
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be interesting to see if the properties of infall changes for lower 
mass (Mhaio < 5 ■ 10^" M©) together with the intrinsic proper- 
ties of galaxies. In addition a systematic study of biases induced by 
the estimators of angular correlations should be conducted, e.g. the 
mass weighted errors we introduced in section l3!2l 

Observationaly, the synthetic halo described in section 15.11 
could be compared to stacked satellite distributions relying on 
galactic surveys such as the SDSS. Once the anisotropy has been 
propagated to the inner regions of the galactic halo following the 
method sketched in section |6| we should be in a position to com- 
pile a synthetic edge-on galactic disk and compare the flaring of 
the disk with the corresponding predictions. The residual preferred 
orientation of galactic disks around more massive objects discussed 
in section l5!2l should be observed on the scales ^ 500kpc. 

Using larger simulations will allow us to combine high resolu- 
tion with the statistics required to detect the anisotropic accretion of 
mass and angular momentum. A wide range of halo masses will be- 
come accessible and the halo mass dependency of our findings will 
be constrained without suffering from the lack of statistics. Better 
angle determinations will naturally follow from a better resolution 
and will improve the accuracy of our quantitative results. Resimu- 
lations (zoom simulations) should give access to a larger range of 
satellite masses, while we were here mostly sensitive to the biggest 
substructures. Large infalling objects are likely to feel differently 
the effects of tidal forces or dynamical friction than smaller satel- 
lites. Resimulated haloes allow us to investigate the dependency on 
the spatial distribution of satellites with their masses correspond- 
ing to a given cosmological environment. However using only a 
few resimulations may not be sufficient to overcome cosmic vari- 
ance and, given the difficulty to produce a large number of high 
resolution haloes, such a project remains challenging. 

The inclusion of gas physics in these simulations and their 
impact on the results is the natural following step. For example, 
gas filaments are known to be narrower than dark matter filaments, 
thus we expect to see a higher level of anisotropy in the distribu- 
tion of accreted gas by the haloes. Furthermore, the transmission 
of angular momentum from one parcel of gas to another (or to the 
underlying dark matter) may be highly effective and would lead 
to higher homogeneity of the properties of the accreted angular 
momentum direction, enhancing the effect of spin alignments. The 
loss of angular momentum from the gaz to the halo will lead to a 
modification of our pure dark matter findings. Yet, the inclusion of 
gas physics in simulations would force us to address issues such as 
the over-cooling, the requirement to take star formation and related 
feedback processes into account. It remains that in the longer term, 
the inclusion of gas physics cannot be avoided and will give new 
insights into the anisotropic accretion of matter by haloes. 
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APPENDIX A: LINEAR RESPONSE OF A SPHERICAL HALO TO INFALLING DARK MATTER FLUXES 

In the followi ng section, we extend to open spherical stellar systems the formalism developed by IXremaine & Weinber j il984) and e.g. 
iMuralil ^99^ by adding a source term to the collisionless Boltzmann equation.^ For an open system, the dark matter dynamics within the 
R200 sphere is governed by the collisionless Boltzmann equation coupled with the Poisson equation: 

^ = ^ + {F,H} = s''{r,^^,t), and V^* = 47rG / d^i;F(v). (Al) 

where { } is the standard Poisson bracket, F(r, v, t) is the system's distribution function submitted to ^(r, t), the total gravitational potential 
(self-gravity + external perturbation). The r.h.s. of JAU is non-zero because of infalling fluxes from the environment which require adding a 
source term, s'^(r, v, t), to the Vlazov equation. We may now discriminate between a stationary part corresponding to the unperturbed state 
from a weak time-dependent perturbation induced by the environnement. Thus the DF can be written as = Fo + /. Provided the mass 
of the incoming flux of dark matter is small compared to the mass of the halo, we may assume that / is small compared to Fq. Similarly, 
the Hamiltonian H of the system can be expanded as Hq + AH, with AH — + ip where tp'^ and tp stand respectively for the external 
perturbative potential and for the small response in potential of the open system. 



Al The Boltzmann equation in action-angle 

Given the periodicity of the system, the most adequate representation of a spherical halo corresponds to action-angle variables iGoldsteii] 
il950h ). The linearized Boltzmann equation in such a representation is: 

(A2) 



^^^^ + »k ■ u^MI, t) = zk . ^ Affk(I, t) + sUl, t). 

The new variables are the actions I and the angles together with the angular rates u! = dw/dt. In eQuation lA2l we have Fourier expanded 
the linearized eguation lAl l over the periodic angles: 



X(r,v,t) = ^Xk(I,t)exp(ik-w), with Xi,{I,t) = ^ dV exp(-ik • w)X(r, v, t), 



(A3) 



where X is any function of (r, v, t) with k being the Fourier triple index corresponding to the three degrees of freedom on the sphere. The 
equilibrium state Fo does not depend on time nor angles since it is assumed to be stationary. Then the solution to <A2t can be written as: 



MI,t)= I drexp(zk.^(r-t)) jk ■ ^ [7/,k(I, r) + ^£(1, r)] + .s^(I, r) 



(A4) 



where we have written A//k(I, t) = V'k(I, t) + 'i/'k(I, ''')• We can integrate <A4t over velocities and sum over k to recover the density 
perturbation: 



p(r, t) = / dr d^v ( exp(ik ■ tj(r — + ik ■ w) 
1. J —00 J 



«k-^[Vk(I,r)+Vk(I,T)] + 4(I,r) 



Let us expand the potential and the density over a bi-orthogonal complete basis function {i/"'"' , p'^'jsuch that 
^{r,t) = ^a„(t)V'["'(r) ; p{r,t) = ^ a„(t)p["l (r) ; VV'"' = 47rGpl"' ; f d^rV-'"'* (r)pW (r) = 

n n 

The external potential can be expanded along the same basis as: 
^^(r,t)=^fo„(t)^["l(r). 



(A5) 



(A6) 



(A7) 



Note that in Eq. <A6t the expansion runs over a triple index n = (n, £, m) corresponding to the radial, azimuthal and alt-azimuthal degrees 
of freedom, while in Eq. <A6> the three coefficients are not independent since the radial variation of the external potential is fixed by its 
boundary value on the sphere i?200. Making use of the biorthogonality, multiplying <A5> by ^^^^*{r) for some given p and integrating over 
r yields: 



1, — oc 



dr / / d^vd^rexp(ik ■ a;(r - t) + zk- w)?/)'P'*(r) 



E zk . ^ [a„(r) + b^ir)] ^"1(1) + 4(1, r) 



(A8) 



This is formally equivalent to summing the response of the halo to a point-like particle for all entering particles. 
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A2 Self-consistency of the response 

We may now swap from position-velocity to angle-action variables since d^vd^r = d^wd^I. In 1A8} only i/''"' (r) depends on w so we may 



carry the w integration over t/j'p'*, yielding which leads to 

(A9) 



k 



CpW = / dr / d^Iexp(jk ■ aj(r - t)) 



. ^ [a„(r) + 6n(r)] + 4(1, r) 



Note that the last term of Eq. <A9> corresponds to the modulated potential along the unperturbed trajectories weighted by the number of 
entering particles with (v, f2) at time r. This is expected since it just reflects the fact that could have linearly summed over all incoming indi- 
vidual particles (since the interaction between particles in a collisionless fluid is purely gravitationnal). In this sense, this terms corresponds 
to a ray tracing problem in a variable index medium. Note also that Eq. <A9t does not account for dynamical friction since we integrate over 
the unperturbed trajectories. At this point, we expand the source term over a complete basis; this basis should also describe (known) velocity 
space variations. We assume that such a basis (^'"' (r, v) exists. We write: 

s"(r,v,t) = Vc„(f)<^'"i(r,v) so 4(1, t) = V c„(r)4"'"(I) where 4"'"(I) = TTr^ / dV exp(-ik • w)<^'"i (r, v). (AlO) 



(27r)a 

Calling a(T) = [ai (r), • • ■ , an(T)], b(r) = [bi(r), 61,(1")], c(r) = [ci (r), • • • , Cn(T)], and 0(r) the Heaviside function, we define 
two tensors: 

ifpn(r) = (1 - e(r)) . ^ /"d^Iexp(zk.a;r)2k.^Vlf'*(I)V'L"'(I), (AH) 
k 

which depends only on the halo equilibrium state via Fo and 

H^t) = (1 - e(r)) . /d'lexp(*k ■ u;r)a|r'^(I) Vlr'*(I), (A12) 
k 

which depends only on the expansion basis, equation <A9> becomes: 

/"OC 

&{t)= I dr{K(r-f)-[a(r) + b(r)] + H(T-t)-c(r)}. (A13) 



We now perform a Fourier transform with respect to time, hence convolutions become multiplications and we get: 

a(p) = (I - K(p))-^ • [K(p) ■ b(p) + H(p) ■ c(p)] , (A14) 
where p stands for the frequency conjugate to time. The computation of the variance-covariance matrix is straightforward: 

(a - a*^) = ((I-K)"^ • [K- b + H-c] • [K- b + H-c]^* ■ (I - K)"^*^) , (A15) 

where I the identity matrix. Note that (a • a* ) involves autocorrelation like (b ■ b* ) and (c ■ c* ) but also cross correlation terms such 
as (b ■ c* ). In other words, recalling that b and c stand respectively for the expansion coeficients of the external potential, Eq. <A7t . and 
the parametrized velocity distribution, Eq. <A10> . their cross-correlation will in fine modify the correlation of the response of the inner halo. 
Two-points statistics are sufficient to caracterize stationnary perturbations and therefore the induced response. Nevertheless, higher statistics 
of the response can be easily expressed in terms of higher order correlations of the pertubation if needed. For example, it can be shown 
that the three-point correlation function of the response's coefficients can be written as function of the two- and three-points correlation of 
the perturbations' coefficients. There are yet quite a few caveats involved; for instance, it is not completely clear today that we have a good 
understanding of what the unperturbed distribution function of a halo-l-disk should be. 



A3 The source term 



A possible choice^ for the source term consistent with the first two velocity moments of the entering matter, involves constructing s"^ (r , v, t) 
in the following manner: 



s'=(r,v,t) = ^Ym(n 



det| {t)/'djp ,m 

it)\ 

= ^ym(n)fo(r-i?200)Crrx(v,t), 



exp 



>(i) 



'OJp,m.{t) 



(A16) 



where m stands for the two harmonic number, {I, m) and Y^{Q,) = Y^iQ,). Here the Dirac function S-o(r — R200) is introduced since we 
measure the source terms at the virial radius. The global form is Gaussian and is constructed using rop,in, ujp^^m, TUpaiaj ,m, the harmonic 



An alternative choice is made in lAubert & PichonI |2004 to account for the bimodality of the velocity distribution. 
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components of respectively the mass flux density field, velocity flux density vector field and the specific kinetic energy flux density tensor 
field measured on the -R200 sphere. When taking the successive moment of this flux distribution over velocity, we get: 

d^vs''(r, v) = TXJp{r) , J d^vvs^(r,v) = ropv(r) , (A17) 

while 

f dV(^, - ^) fa- - ^)s-(r,v) = ^p...,(r) + ( Vy^(n)^(r ~ i?2oo) '^f'-'if - ) « (r), (A18) 

J lA^p ZUp '(^p,m{t) i^pir) / 

so that the Ansatz, Eq. <A16t , satisfy the first two moments, and approximatively the third moment of the fluid equations. Let us now expand 
Cm(v, t) over a linear complete basis, say b-splines covering the radial velocity component and spherical harmonics for the angle distribution 
of the velocity vector: 

Cm(v,t) = ^C^,c«W&a(v). (A19) 

The particular choice of Eq. <A16> has led to the parametrisation: 

Cn(t) = C^,a(t) and 0'"'(r, v) = bc,(w)Y^{n)5u{r - i?20o) , (A20) 
while Eq. <A10> becomes 
1 

(2V)3 

Note that we can make use of the Su function occurring in Eq. jA21> since Wr = Wr{T, I). Therefore Eq. jA21t reads: 

'"'■exp(-lk- W) Ym[ft{I,W, l5r[i?200,I])lba(v[I,W, Wr{R200,i)]) , ■/ ^''^^'^yM '^^P i^l-^r ' iSr [-R2OO , I] ) (A22) 



= j—r^ I d'wexp(-ik- w)y„[n(I,w)]6c(v[I,w])5D(r-(I,w)-7?20o). (A21) 



(27r)^ |r(iJ:2oo,I) 

In Eq. <A22t we sum over all intersections of the orbit I with the R200 sphere, at the radial phase corresponding to that intersection (with a 
weight corresponding to ujr/\r\). 

Given Eq. <A6t . Eq. <A16> . together with Eq. <A22t . Eq. <A13t can be recast formally as 

p(r,t) = R{Fo,t,r,n} [tP^{n,T),zup(n,T),zup^{n,T),vjp^^^^{n,T)] (A23) 

which corresponds to the form given in the main text in Eguation lSlI It should be emphasized once again that the splitting of the gravitationnal 
field into two components, one outside of i?200, and one inside, via point particles obeying the distribution Se(r, v, t) is completely arbitrary 
from the point of view of the dynamics. In fact, one should account that (il , i ) should be switched on long before any particles enter i?200 
since no particle is created at the boundary. This last constraint is clearly satisfied by our simulations. 
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APPENDIX B: ADAPTAHOP: A SUBSTRUCTURE 
FINDER BASED ON SADDLE POINT HANDLING 

Dark matter haloes can contain a hierarchy of subhaloes, which 
can be viewed as a tree of structures and sub-structures. Given a 
mass resolution (a finite number of particles such as in our TV-body 
simulations), there is a limit to this hierarchy, which can be for- 
malised as an ensemble of leaves in a tree. The goal here is to draw 
this tree b y applying the simplest principles of Morse theory (e.g. 
|joslH2od2 ')). Morse theory basically involves relating the topology 
of an excursion, e.g., the regions of space with density above a 
given threshold, p > pt , to the set of critical points it contains, 
{x, Vp(x) = 0}, and to the field lines connecting these points 
together, i.e. the curves obtained by following the gradient of the 
density field. In that approach, the smallest substructures, which 
are the leaves of the tree, can be identified as peak patches, i.e. en- 
sembles of field lines converging to the same local maximum. The 
connectivity between substructures is ruled by the saddle points, 
which are local maxima in the surfaces defining the contours of the 
peak patches: from the knowledge of these saddle points and the 
local maxima they connect, it is possible to extract the full tree of 
structures (haloes) and sub-structures (subhaloes) in four steps: 

(i) In order to eliminate, at least partly, the effects of Poisson 
noise and to have an estimate of the local density as close as pos- 
sible to a Morse function,'^ while conserving as much as possible 
details of the distribution, we perform adaptive smoothing of this 
distribution wit h sta ndard SPH technique (Smooth Particle Hydro- 
dynamics, e.g. iMon aghan (1992)). This smoothing assumes that 
each particle is a smooth spherical cloud of given radius R, e.g. a 
spline Sir). For each particle, the list of its A^sph closest neighbors 
is found, typically TVsph of a few tens (here we take A^sph = 64). 
The distance from the furthest neighbor fixes R, while the SPH 
density at the particle of interest is estimated by a summation over 
its neighbors with weight S(r). To find rapidly the closest neigh- 
bors of each particle, we use a standard Oct-tree algorithm, which 
decomposes hierarchically space in subcells untill they contain zero 
or one particle. 

(ii) The leaves of the tree of structures and substructures are 
identified while associating each particle to the peak patch it be- 
longs to. This is performed by a simple walk from particle to parti- 
cle, while following the gradient until convergence: at each step of 
the walk, the SPH density of the particle is compared to its A^hop 
closest neighbors (which were stored during the SPH smoothing 
step), the particle for the next step of the walk being the one with 
the largest SPH density . We take A'^hop = 16, as advocated by 
lEisenstein&HulU 19981) . 

(iii) For each leave of the tree, the connections with the other 
leaves are created by searching the saddle points on the intersect- 
ing surfaces Sij between peak patches i and j. Each surface Sij is 
made of particles belonging to one of the peak patches and having 
at least one of their closest neighbors among A'hop in the other 
peak patch, and vice versa. If the set Sij contains only particles 
belonging to i or only particles belonging to j, the connection be- 
tween i and j is considered as non significant (because non sym- 
metric) and eliminated, Sij = . Saddle points are local maxima in 
Sij. To establish the connectivity as a function of a density thresh- 
old, only the highest saddle point matters, when there are several. 



^ i.e. a smooth function such that the ensemble of critical point is discrete 
and the matrix of second derivatives in their neighborhood is non degener- 
ate. 



The search for this saddle point involves finding the maximum of 
the SPH density among particles belonging to Sij. We proceed as 
follows to estimate accurately the SPH density in Sij . For each par- 
ticle Am Sij, say belonging to peak patch i and with density pA, 
we consider the list of its closest neigbors among A^hop belong- 



ing to peak patch j, with density p^, k 



The density associated to this particle in Sij is then given by 
p = min(p,4,pfc). By applying this procedure, we locate accu- 
rately Sij and avoid slight overestimation of the SPH density at the 
saddle point. 

(iv) It is possible to build the tree of structures and sub- 
structures, when the list of neighboring leaves to which a given 
leave is connected is given, as well as the corresponding saddle 
points. This is performed recursively by increasing progressively 
a threshold parameter, pt, from an initial value, pth, correspond- 
ing to the typical overdensity used to select galaxy haloes, here 
called structures. A typical choice for pth is Pth ~ 81, which 
corresponds approximately to frie nd-of -friend haloes selected with 
a linking paremeter b = 0.2 (e.g.. lEisenstein & Hu Ml998.) ). Sup- 
pose we are at step n of the process and let us compute step n + 1. 
At this point, we are sitting on a branch of the tree -a structure 
or a sub-structure- and we aim to draw the details of this branch. 
This (sub-) structure contains a number of peak patches connected 
by saddle points of densities ps. For the considered value of pt, the 
connections inside that (sub-) structure are examined and destroyed 
when Ps < Pt. The (sub-) structure is then broken into as many 
components as necessary. During the process, the particles above 
Pt belonging to each sub-component are tagged, which allows us 
to determine at any time various properties of a given (sub-) struc- 
ture, namely the number of particles it contains, its mass, its aver- 
age and maximum SPH density, for possible application to various 
morphological criterions of selection. One such criterion is Poisson 
noise. In order to asses if a given sub-structure containing A'^ par- 
ticles should be considered as statistically significant compared to 
Poisson noise, its average density must be sufficiently significant 
compared to pt : 



(p)sub — structure ^ Pt 



1 + 



/Pois 



(Bl) 



where /poisson is a "/poissoncr" detection parameter, typically a few 
unities. A good choice is /poisson = 4. If the sub-structure is below 
this threshold, it disappears, i.e. it is not considered in next step of 
the recursion. At the end of the selection, two situations are possi- 
ble: (i) two sub-structures or more are detected and new nodes are 
created in the tree (ii) the (sub)-structure was not broken into mul- 
tiple components and nothing happens at this step. The process is 
then repeated on the new sub-structures by increasing locally the 
threshold pt : 



pt 



Pt X 



1 



(B2) 



until there is only one peak patch in the (sub-)structure. Note that 
the Poisson noise selection, eq. <B1> is not applied to the haloes, 
when Pt = PTH. 

At the end of the process, one obtains a tree of structures and 
sub-structures, where each node of the tree corresponding to a (sub- 
)structure, with its position, its number of particles, is mean square 
radius, its average and maximum SPH density, and the density ps of 
the highest saddle point which connects it to another substructure. 
In addition, a flag is given to each particle. This flag is a pointer to 
the node the closest possible to a leave (if not a leave), which allows 
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one to find recursively the list of particles belonging to any (sub- 
)structure and thus perform some more ellaborate post-treatement, 
such as some relying on dynamical prescriptions (boundness). The 
difficulty in that case is to estimate accuratly the gravitational po- 
tential. Its computation can be rather costly, since "pealing" the 
(sub-)haloes requires iterating several forward and backward walks 
in the tree of structures and sub-structures with corresponding cal- 
culations of the gravitational potential. Our prescription is there- 
fore at the present time purely morphological and does not involve 
the estimate of the gravitational potential. The current implemen- 
tation is rather fast, most of the CPU time being taken by the SPH 
smoothing, e.g. 1-3 hours on 16 millions particles on current fast 
scalar processors. 

Our algo rithm is called ADAPTA HOP since we aim at im- 
proving HOP iEisenstein & Hull il998t) ): the first two step above 
are exactly the same as in HOP, but the last two are different. In- 
deed, in HOP, the idea is to combine informations on the saddle 
point densities, ps, on the local maxima, /Omax, inside a connected 
set of peak patches to decide whether it has to be broken into mul- 
tiple disjoint haloes. The aim of HOP is indeed to improve stan- 
dard friend-of-friend methods in order to obtain more compact and 
spherical haloes. The goal of ADAPTAHOP is quite different since 
it focuses on sub-structure detection. 

In spirit, ADAPTAHOP is in fact very similar to the subst ruc- 
ture fi nder of lSpringell (l99i): SUBFIND (see also Sorinael et^d] 
i200ll) '). Of course, there is a major difference, since SUBFIND 
has in addition a sophisticated dynamical prescription involving 
exact calculation of the gravitationnal potential. Springel uses also 
a slightly more elegant method to construct the tree of structures 
and substructures prior to dynamical post-treatment. After step one 
above, the idea is to rank the particles by decreasing density and 
treat them in this order. Investigating the distribution of particules 
is such a way is equivalent to examining isocontours of decreasing 
density. It uses (as in ADAPTAHOP) the closests neighbors of a 
particle to decide if the particle examined during the process (i) cre- 
ates a new (sub-)structure since it is isolated (ii) belongs to an ex- 
isting substructure or (iii) connects two substructures, which makes 
the construction of the tree of structures and substructures much 
simpler than in ADAPTAHOP and more accurate, since there no 
need for using the threshold parameter pt. In SUBFIND, no treat- 
ment is made to account for the local Poisson noise: it is not nec- 
essary because of the dynamical post-processing, which destroys 
unbounded structures. 

It is important to note that since ADAPTAHOP has no dynam- 
ical post-treatment, it gives slightly different results compared to 
SUBFIND in its present form. In partic ular, for a given suffici ently 
massive dark matter halo , SUBFIND jSpringel etlH j200ll) ') de- 
scribes it in terms of a large, smooth central component, and a 
bunch of much less massive sub-haloes. In ADAPTAHOP, the re- 
sult is quite similar, except that the central component is much less 
spatially extended (it is extended up to the isocontour level corre- 
sponding to the saddle point connecting it to a sub-halo), and it is 
therefore less massive. 

Fig. lBll illustrates how well ADAPTAHOP performs in one of 
the simulations we realized for this work, for the most massive halo 
detected in this realization. 



APPENDIX C: STATISTICS ON THE SPHERE 

When dealing with spherical fields, there are different ways to char- 
acterize their angular structure. In the present paper, we essentially 




Figure Bl. Illustration of the output of ADAPTAHOP for one of the simu- 
lations of this work. A sphere of radius 5 Mpc centered on the most massive 
halo is represented. In the upper pannel, the dark matter density is shown 
using a logarithmic scale. Darker regions correspond to higher density con- 
trasts. The lower pannel displays the detected subhaloes (i.e. the most el- 
ementary structures con'esponding to the peak patches or the leaves of the 
tree). The size of the circle scales with M^/^, where M is the mass of the 
subhalo. Most of the subhaloes seen on the figure belong to the most mas- 
sive halo. Clearly, ADAPTAHOP is rather successful! at detecting all the 
significant substructures. 



deal with centered statistics, i.e. we describe the angular structura- 
tion of scalar or vector fields relative to a specific direction, defined 
by the halo or satellite's spin S. Let us first formally introduce fil- 
tering on the sphere, statistical and angular averages, and present 
one-point statistics (probability distribution functions) while post- 
poning two-point statis tics (correlation function s, or excess proba- 
bility of joint events) to lAubert & PichonI i2004) . 
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Cl One point statistics 

For any field, x, on the sphere, let us introduce the smoothed field, 
{x)a (filtered on scale a) , as 



Note that we can think of the harmonic coefficients, a™ in- 
troduced in section I7!4l as a specific type of filtering, where the 
window function, Wa, is replaced by an axisymmetric spherical 
harmonic, Y"{fl): 



(X)c(n) 



—, /ec(fi'-fi)x(fi')dJ^' (Cl) , 1 /"^o^/^x ,^ a°e 

VJ J ■CO 



(C2) 



where Bq stands for the top hat function, 

e<,(n) = 1 if \-d\ < a, (C3) 

and Wa is defined bv lC7l the standard top hat filter on the sphere. 

Consider now the centered top-hat-filtered (on scale a) field, 
[x\a, defined by 



[x\a = {x)c{-), 



(C4) 



\^7^ j e„(^'-V2)x(n')dfl', (C5) 



Wa[n')x{n')An' 



(C6) 



Note that Eq. <C6> defines Wa- Our filtering is now centred, in that 
the average is carried on a window which is centred at the equatorial 
plane (since we are in this paper interested in the polarisation of 
accretion processes with respect to that plane). Let us also introduce 
the average of x on the sphere, as 



X = (x)^/2 = 



x{n)An, 



We may also for a given x define its contrast as 



(C8) 



S^ = --l. 

X 

Note that, in contrast to standard cosmology, we expect that x 7^ 
{x), (i.e. no ergodicity) since the angular average over one virial 
sphere is not representative of the whole cosmological set, and 
since {x) depends on -d whereas x doesn't. As a consequence, 

{5.) = (-) - 1 / - 1 . 
X [x) 

Consider now the top-hat-filtered-centred flux density contrast, 
defined by 

[Su,]a = {S^)a(n/2) = ^ J Wa^VJ - 1 . (C9) 

Since by construction, is a filtered version of 5-^, it inherits 

some of it statistical properties. In particular, the PDF of 5t^{tt/2) 
and [S^]a should be quite similar provided a is small enough. 

In the main text, we consider the anisotropic parameter, Sm = 
[Spv^]„/s, which therefore corresponds formally to the centered 
top-hat-smoothed (on scales of vr/S) mass flux density contrast. 
Following the same spirit, we could also consider quantities such 
as [Sp^^y2]j^/s, which would measure the anisotropy in the accreted 
kinetic energy: the excess of accreted kinetic energy should allow 
us to track the excess of incoming virialised objects in the equa- 
torial plane without performing their explicit identification. One 
should also consider [Spv,.l]tt/s, the anisotropy in the accreted mo- 
mentum, since this quantity is directly related to the torque applied 
onto the system by the infall. More generally still we could inves- 
tigate {Sui)a{'&), the flux density contrast top-hat-smoothed on a 
ring of size a centred on i?. 



We can also write zu in terms of spherical harmonics 
1 /■ ,„ 1 



zu = — / zudfl = — r= 
47r I 



47r 



Therefore we obtain : 



sign(a[j)^C^ 



(CIO) 



(CU) 



(C12) 



where Co = !aoiV4 

TT is the ^ = component of the angular power 

spectrum Ci. 

Since a step function can be expanded along spherical har- 
monics as 



therefore, [5tc]a defined by Ea. lC9l obevs 
Taking x = pv^ for example, we have : 



(C13) 



(C14) 



(C15) 



where 



piv = — /di9d<^pu,.(i9,V9)sini? = (C16) 
47!- } V47r 

Since 

y di9d(^y/"(i9, (/p) sin-i? = yT^&m&rai^ 

(e.g. lVarshalovich et alJ il983) ). we find : 

= al - %/47r5!o^mO- (C17) 
We finally obtain : 

= ^aryr(^,^) - 1. (C18) 

APPENDIX D: CONVERGENCE ISSUES 
DI Substructures & the haloe's spin 

For each tree of substructure-satellites, we computed the total spin 
inside the mother structure, '&m, and the momentum of each sub- 
structure inside the mother structure, Ls. Then we compared the in- 
ner satellites' and the contribution of the core to the mother's spin. 
The comparison is only made on the components of the substruc- 
tures momentum parallel to Sm- The results are shown in Fig. lDll 
We plotted the total contribution of satellites to the mother's spin 
versus the core's contribution. From the barycenter of the distribu- 
tion shown in Fig. lDll it appears that substructures contain about 
80 % of the total host's spin with a satellites contribution of 50% 
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Figure D2. Comparison of {&m) for different classes of halo mass at z=0. 
The error bars stand for the 3(t error. The thin lines separate the three classes 
of mass: 5 ■ lO^^A/g < m < 1.25 ■ 10^^ Mq, 1.25 ■ 10^^ Mq < m < 
2.5 ■ IO^^Mq and m > 2.5 ■ lO^^M©. Each class contains 16500 haloes. 



Consequently, no class of mass dominates when all the haloes are 
being used in the computation of {5m). 



Figure Dl. Comparison of the substructure's and the core's contribution to 
the amplitude of the mother's spin and to the mother's mass. Top: Compar- 
ison of the core's contribution to the mother's spin compared to the con- 
tribution of all the satellites for each mother detected in our simulations. 
Bottom: same comparison but for the core's and satellites' mass relative to 
the mother's total mass. In both figures, the symbol indicates the barycenter 
of the cloud of points while the thick line's slope is unity. While the total 
mass is dominated by the core's contribution, the mother's spin is domi- 
nated by satellites showing that their specific orbital momentum is more 
important than that of the core. 



and about 30 % for the core. The bottom panel shows the total con- 
tribution of substructures to the mother's mass versus the contribu- 
tion of the core. As expected given the definition of the core, we 
found that the relative proportion are almost reversed compared to 
the previous plot. A core contains about half of the total mass while 
satellites represent about 40% of the total mass. Clearly the specific 
angular momentum is larger in satellites than in the core. The dis- 
tance of satellites relative to the mother's centre and their velocities 
induce a 'lever arm' effect. Even if satellite remnants are light in 
terms of mass they are important if not dominant for the spin of the 
galactic system. This effect also suggests that the mother's spin is 
aligned with the orbital momentum of infalling satellites because 
they determine the direction of the halo's spin. 



D2 the mass dependence of {5m) 

We measured the average excess of accretion {5m) ( see sec- 
tion |33j for three different class of masses at redshift z=0: 5 • 
lO^^M© < m < 1.25 ■ 1O"M0, 1.25 • IO^^Mq < m < 
2.5 ■ IO^^Mq and m > 2.5 • lO^^M©. Each class contains ap- 
proximatively 16 500 haloes. The results are shown in Fig. lD2l It 
is found that {5m) increases with mass but does not change signif- 
icantly even if the three classes cover different mass magnitudes. 



